home *** CD-ROM | disk | FTP | other *** search
/ Night Owl 6 / Night Owl's Shareware - PDSI-006 - Night Owl Corp (1990).iso / 006a / acm525.zip / ACM525.FOR
Text File  |  1991-02-15  |  90KB  |  2,256 lines
  1. C     ALGORITHM 525, COLLECTED ALGORITHMS FROM ACM. THIS WORK
  2. C     PUBLISHED IN TRANS. MATH. SOFTWARE, 4(1), PP.82-94.
  3.  
  4. C
  5. C  **** THIS PROGRAM RUNS ADAPT WITH DATA READ IN FROM CARDS
  6. C       IT IS DESIGNED TO TEST ADAPT WITH THE 20 FUNCTIONS IN F(X).
  7. C       IT READS DATA FOR EVERYTHING EXCEPT EDIST(=ATYPE) AND PRINTS
  8. C       A HEADING WITH THE FUNCTIONS NAME. IT CALLS AND TIMES ADAPT,
  9. C       THEN INDEPENDENTLY CHECKS THE ERROR IN ADAPT AND PLOTS THE
  10. C       FUNCTION AND ERROR CURVE.
  11. C
  12. C NOTES -- USER MUST SUPPLY SYSTEM ROUTINES SECOND AND GRAPH
  13. C          INTEGRATION USES ERRINT FROM ADAPT.
  14. C          PRESENT DIMENSIONS (ON ZKNOTS + ZCOEFS) PREVENT USING
  15. C             ERRINT TO CHECK ACCURACY FOR MORE THAN 100 KNOTS OR
  16. C             POLYNOMIAL DEGREE 6.
  17. C          QPOLY USED BY ERRINT IS A SINGLE ARGUMENT VERSION OF PPOLY.
  18. C          SEE COMMENTS IN ADAPT ABOUT DIMENSION CONSTRAINTS AND
  19. C             PORTABILITY.
  20. C          IF ADAPT AND THIS DRIVER ARE CHANGED TO DOUBLE PRECISION
  21. C             ONE PROBABLY WANTS TO LEAVE TIME,TSTART,TSTOP,XVALU,
  22. C             GRAF1,GRAF2 AS SINGLE PRECISION
  23. C             THEY ARE SEPARATELY DECLARED HERE.
  24. C
  25. C
  26.       REAL A, ACCUR, B, BLEFT, BRIGHT, CHARF, DDTEMP, DSCTOL, ERROR,
  27.      * ERRORI, FACTOR, FINTRP, FLEFT, FRIGHT, NORM, XBREAK, XDD,
  28.      * XINTRP, XLEFT, XRIGHT
  29.       DIMENSION XBREAK(20), DBREAK(20), BLEFT(20), BRIGHT(20)
  30.       DIMENSION XLEFT(50), XRIGHT(50), FACTOR(12), FMESGE(6)
  31.       DIMENSION DDTEMP(12,12), FINTRP(10), FLEFT(6), FRIGHT(6),
  32.      * XDD(12), XINTRP(10)
  33.       INTEGER BOTH, BREAK, DBREAK, DEGREE, EDIST, FMESGE, RIGHT,
  34.      * RIGHTX, SMOOTH
  35.       LOGICAL DISCRD, FATAL, FINISH
  36.       DIMENSION XKNOTS(140), COEFS(140,12)
  37.       REAL COEFS, ECHECK, F, XKNOTS
  38.       REAL TIME, TSTART
  39.       DIMENSION NAMES(25,10)
  40.       EXTERNAL F
  41.       COMMON /INPUTZ/ A, B, ACCUR, NORM, CHARF, XBREAK, BLEFT, BRIGHT,
  42.      * DBREAK, DEGREE, SMOOTH, LEVEL, EDIST, NBREAK, KNTDIM, NPARDM
  43.       COMMON /RESULZ/ ERROR, KNOTS
  44.       COMMON /KONTRL/ DSCTOL, ERRORI, XLEFT, XRIGHT, BREAK, BOTH,
  45.      * FACTOR, FMESGE, IBREAK, INTERP, LEFT, MAXAUX, MAXKNT, MAXPAR,
  46.      * MAXSTK, NPAR, NSTACK, RIGHT, DISCRD, FATAL, FINISH
  47.       COMMON /COMDIF/ DDTEMP, FINTRP, FLEFT, FRIGHT, XDD, XINTRP,
  48.      * LEFTX, NINTRP, RIGHTX
  49.       COMMON /SAVEIT/ IKNOT
  50. C
  51.       COMMON /FDATA/ JFUNK, KOUNT, MAXK
  52.       COMMON /INFO/ ECHECK, KFUNK, TIME, TSTART
  53.       DATA NAMES(1,1), NAMES(1,2), NAMES(1,3), NAMES(1,4), NAMES(1,5),
  54.      * NAMES(1,6), NAMES(1,7), NAMES(1,8), NAMES(1,9), NAMES(1,10) /
  55.      * 4HSIMP,4HLE A,4HLGEB,4HRAIC,4H SIN,4HGULA,4HRITY,4H AT ,4H0.0 ,
  56.      * 4H    /
  57.       DATA NAMES(2,1), NAMES(2,2), NAMES(2,3), NAMES(2,4), NAMES(2,5),
  58.      * NAMES(2,6), NAMES(2,7), NAMES(2,8), NAMES(2,9), NAMES(2,10) /
  59.      * 4HINTE,4HRIOR,4H SIN,4HGULA,4HRITY,4H AT ,4HX = ,4H.307,4H7070,
  60.      * 4H7707/
  61.       DATA NAMES(3,1), NAMES(3,2), NAMES(3,3), NAMES(3,4), NAMES(3,5),
  62.      * NAMES(3,6), NAMES(3,7), NAMES(3,8), NAMES(3,9), NAMES(3,10) /
  63.      * 4HBROK,4HEN L,4HINE ,4H- BR,4HEAKS,4H AT ,4H1,2,,4H3.5 ,4HAND ,
  64.      * 4H4   /
  65.       DATA NAMES(4,1), NAMES(4,2), NAMES(4,3), NAMES(4,4), NAMES(4,5),
  66.      * NAMES(4,6), NAMES(4,7), NAMES(4,8), NAMES(4,9), NAMES(4,10) /
  67.      * 4HHUMP,4H AT ,4H2.5 ,4HRECI,4HPROC,4HAL O,4HF QU,4HARTI,4HC   ,
  68.      * 4H    /
  69.       DATA NAMES(5,1), NAMES(5,2), NAMES(5,3), NAMES(5,4), NAMES(5,5),
  70.      * NAMES(5,6), NAMES(5,7), NAMES(5,8), NAMES(5,9), NAMES(5,10) /
  71.      * 4HSIMP,4HLE C,4HUBIC,4H = (,4H4X-2,4H)**3,4H - (,4H4X-2,4H)   ,
  72.      * 4H    /
  73.       DATA NAMES(6,1), NAMES(6,2), NAMES(6,3), NAMES(6,4), NAMES(6,5),
  74.      * NAMES(6,6), NAMES(6,7), NAMES(6,8), NAMES(6,9), NAMES(6,10) /
  75.      * 4HVARI,4HABLE,4H OSC,4HILLA,4HTION,4H FRO,4HM SI,4HN(X*,4H*2-1,
  76.      * 4H)   /
  77.       DATA NAMES(7,1), NAMES(7,2), NAMES(7,3), NAMES(7,4), NAMES(7,5),
  78.      * NAMES(7,6), NAMES(7,7), NAMES(7,8), NAMES(7,9), NAMES(7,10) /
  79.      * 4HEXPO,4HNENT,4HIAL ,4H    ,4H    ,4H    ,4H    ,4H    ,4H    ,
  80.      * 4H    /
  81.       DATA NAMES(8,1), NAMES(8,2), NAMES(8,3), NAMES(8,4), NAMES(8,5),
  82.      * NAMES(8,6), NAMES(8,7), NAMES(8,8), NAMES(8,9), NAMES(8,10) /
  83.      * 4HFIFT,4HH DE,4HGREE,4H POL,4HY PL,4HUS F,4HINE ,4HSAW ,4HTOOT,
  84.      * 4HH   /
  85.       DATA NAMES(9,1), NAMES(9,2), NAMES(9,3), NAMES(9,4), NAMES(9,5),
  86.      * NAMES(9,6), NAMES(9,7), NAMES(9,8), NAMES(9,9), NAMES(9,10) /
  87.      * 4HSEVE,4HNTH ,4HDEGR,4HEE P,4HOLYN,4HOMIA,4HL   ,4H    ,4H    ,
  88.      * 4H    /
  89.       DATA NAMES(10,1), NAMES(10,2), NAMES(10,3), NAMES(10,4),
  90.      * NAMES(10,5), NAMES(10,6), NAMES(10,7), NAMES(10,8), NAMES(10,9),
  91.      * NAMES(10,10) /4HWILD,4H VAR,4HIATI,4HON B,4HETWE,4HEN .,4H47 A,
  92.      * 4HND .,4H48  ,4H    /
  93.       DATA NAMES(11,1), NAMES(11,2), NAMES(11,3), NAMES(11,4),
  94.      * NAMES(11,5), NAMES(11,6), NAMES(11,7), NAMES(11,8), NAMES(11,9),
  95.      * NAMES(11,10) /4HLOTS,4H OF ,4HOSCI,4HLLAT,4HION ,4H - F,4HIRST,
  96.      * 4H EXA,4HMPLE,4H    /
  97.       DATA NAMES(12,1), NAMES(12,2), NAMES(12,3), NAMES(12,4),
  98.      * NAMES(12,5), NAMES(12,6), NAMES(12,7), NAMES(12,8), NAMES(12,9),
  99.      * NAMES(12,10) /4HLOTS,4H OF ,4HOSCI,4HLLAT,4HION ,4H - S,4HECON,
  100.      * 4HD EX,4HAMPL,4HE   /
  101.       DATA NAMES(13,1), NAMES(13,2), NAMES(13,3), NAMES(13,4),
  102.      * NAMES(13,5), NAMES(13,6), NAMES(13,7), NAMES(13,8), NAMES(13,9),
  103.      * NAMES(13,10) /4HSING,4HULAR,4HITIE,4HS AT,4H BOT,4HH EN,4HD PO,
  104.      * 4HINTS,4H    ,4H    /
  105.       DATA NAMES(14,1), NAMES(14,2), NAMES(14,3), NAMES(14,4),
  106.      * NAMES(14,5), NAMES(14,6), NAMES(14,7), NAMES(14,8), NAMES(14,9),
  107.      * NAMES(14,10) /4HRAND,4HOM N,4HUMBE,4HRS A,4HDDED,4H TO ,4HSIN(,
  108.      * 4HX)  ,4H    ,4H    /
  109.       DATA NAMES(15,1), NAMES(15,2), NAMES(15,3), NAMES(15,4),
  110.      * NAMES(15,5), NAMES(15,6), NAMES(15,7), NAMES(15,8), NAMES(15,9),
  111.      * NAMES(15,10) /4HLARG,4HE FU,4HNCTI,4HON -,4H ADJ,4HUSTA,4HBLE ,
  112.      * 4HSIZE,4H    ,4H    /
  113.       DATA NAMES(16,1), NAMES(16,2), NAMES(16,3), NAMES(16,4),
  114.      * NAMES(16,5), NAMES(16,6), NAMES(16,7), NAMES(16,8), NAMES(16,9),
  115.      * NAMES(16,10) /4HSHIF,4HTED ,4HORIG,4HIN -,4H ADJ,4HUSTA,4HBLE ,
  116.      * 4HSHIF,4HT   ,4H    /
  117.       DATA NAMES(17,1), NAMES(17,2), NAMES(17,3), NAMES(17,4),
  118.      * NAMES(17,5), NAMES(17,6), NAMES(17,7), NAMES(17,8), NAMES(17,9),
  119.      * NAMES(17,10) /4HNUME,4HRICA,4HL DI,4HFFER,4HENTI,4HATIO,4HN OF,
  120.      * 4H F(X,4H)   ,4H    /
  121.       DATA NAMES(18,1), NAMES(18,2), NAMES(18,3), NAMES(18,4),
  122.      * NAMES(18,5), NAMES(18,6), NAMES(18,7), NAMES(18,8), NAMES(18,9),
  123.      * NAMES(18,10) /4HCOMP,4HLICA,4HTED ,4HFUNC,4HTION,4H WIT,4HH 7 ,
  124.      * 4HPIEC,4HES  ,4H    /
  125.       DATA NAMES(19,1), NAMES(19,2), NAMES(19,3), NAMES(19,4),
  126.      * NAMES(19,5), NAMES(19,6), NAMES(19,7), NAMES(19,8), NAMES(19,9),
  127.      * NAMES(19,10) /4HJUMP,4H DIS,4HCONT,4HINUI,4HTY A,4HT VA,4HRIAB,
  128.      * 4HLE P,4HOINT,4H    /
  129.       DATA NAMES(20,1), NAMES(20,2), NAMES(20,3), NAMES(20,4),
  130.      * NAMES(20,5), NAMES(20,6), NAMES(20,7), NAMES(20,8), NAMES(20,9),
  131.      * NAMES(20,10) /4HINFI,4HNITE,4H SIN,4HGULA,4HRITY,4H AT ,4HP10A,
  132.      * 4H    ,4H    ,4H    /
  133. C
  134.    10 READ (5,99999) JFUNK, SMOOTH, DEGREE, LEVEL, NORM, A, B, ACCUR,
  135.      * CHARF, NBREAK
  136.       IF (JFUNK.EQ.0) STOP
  137.       IF (CHARF.EQ.0.) CHARF = (B-A)*1.001
  138.       IF (NBREAK.EQ.0) GO TO 20
  139.       READ (5,99998) (DBREAK(K),XBREAK(K),BLEFT(K),BRIGHT(K),K=1,NBREAK)
  140.    20 CONTINUE
  141.       MAXK = 3200
  142.       EDIST = 0
  143.       WRITE (6,99997) JFUNK, (NAMES(JFUNK,K),K=1,10)
  144.       KOUNT = 0
  145. C
  146. C     *****  CALL ON CLOCK TO INITIALIZE  TSTART  IN SECONDS  *****
  147. C
  148.       TSTART = 0.0
  149.       CALL ADAPT(F, XKNOTS, COEFS, 140, 12)
  150.       CALL SUMRY2(XKNOTS, COEFS, 140, 12)
  151.       GO TO 10
  152. 99999 FORMAT (4I3, F4.3, 2F10.3, E12.4, F10.3, 2I3)
  153. 99998 FORMAT (2(I5, 3F11.5))
  154. 99997 FORMAT (//5(3H+++), 18HADAPT FOR FUNCTION, I3, 3H = , 10A4, 4X,
  155.      * 5(3H+++))
  156.       END
  157.       BLOCK DATA
  158. C         SET PARAMETERS FOR FUNCTIONS IN F
  159.       REAL PAR2, PAR3A, PAR3B, PAR4, PAR4A, P5A, P5B, P6A, P6B, P7,
  160.      * P8A, P8B, P8C, P8D, P8E, P8F, P8G, P8H, P9, P10A, P10B, P10C
  161.       COMMON /FPARS/ PAR2, PAR3A, PAR3B, PAR4, PAR4A, P5A, P5B, P6A,
  162.      * P6B, P7, P8A, P8B, P8C, P8D, P8E, P8F, P8G, P8H, P9, P10A, P10B,
  163.      * P10C
  164.       DATA PAR2, PAR3A, PAR3B, PAR4, PAR4A, P5A, P5B
  165.      * /.02,.4,.6,.001,6.,8000.,2037./
  166.       DATA P6A, P6B, P7, P8A, P8B, P8C, P8D, P8E /4000.,3998.,.0001,-2.,
  167.      * 1.5,7.,5.,9.6/
  168.       DATA P8F, P8G, P8H, P9, P10A, P10B, P10C /.6,14.,16.,1.0,.25,.5,
  169.      * .6667/
  170.       END
  171.       REAL FUNCTION F(X, FDERV)
  172. C
  173. C **  AN ARRAY OF TWENTY TEST FUNCTIONS WITH 2 TO 5 DERIVATIVES.
  174. C     THE DERIVATIVE VALUES ARE PLACED IN FDERV - DIMENSIONED AT 5
  175. C     THE FUNCTIONS OF THE SECOND GROUP ARE PARAMETERIZED IN VARIOUS
  176. C     WAYS.  THE PARAMETERS ARE SET IN BLOCK DATA AND AVAILABLE
  177. C     THROUGH THE COMMON BLOCK / FPARS /
  178. C         REQUIRED INPUT INFORMATION IN COMMON BLOCK / FDATA /
  179. C            JFUNK  = INDEX OF THE FUNCTION SELECTED
  180. C            KOUNT  = NUMBER OF F-EVALUATIONS
  181. C            MAXK   = MAXIMUN ALLOWED VALUE OF KOUNT
  182. C         REQUIRED CONTROL INFORMATION IN COMMON BLOCK / KONTROL /
  183. C            FINISH = SWITCH THAT STOPS ADAPT
  184. C
  185.       REAL A, ACCUR, B, BLEFT, BRIGHT, CHARF, DDTEMP, DSCTOL, ERROR,
  186.      * ERRORI, FACTOR, FINTRP, FLEFT, FRIGHT, NORM, XBREAK, XDD,
  187.      * XINTRP, XLEFT, XRIGHT
  188.       DIMENSION XBREAK(20), DBREAK(20), BLEFT(20), BRIGHT(20)
  189.       DIMENSION XLEFT(50), XRIGHT(50), FACTOR(12), FMESGE(6)
  190.       DIMENSION DDTEMP(12,12), FINTRP(10), FLEFT(6), FRIGHT(6),
  191.      * XDD(12), XINTRP(10)
  192.       INTEGER BOTH, BREAK, DBREAK, DEGREE, EDIST, FMESGE, RIGHT,
  193.      * RIGHTX, SMOOTH
  194.       LOGICAL DISCRD, FATAL, FINISH
  195.       REAL C, CP, CS, C1, DDR, DDT1, DDT2, DDT3, DDT4, DDT5, DDT6,
  196.      * DDT7, DDXEP, DL, DR, DS, DT1, DT2, DT3, DT4, DT5, DT6, DT7,
  197.      * DXEP, D2L, D2R, D2S, D3L, D3R, D3S, EPS, EX, FDERV, FD1, FD2,
  198.      * FK, FL, FLL, FR, FRR, F27, PAR2, PAR3A, PAR3B, PAR4, PAR4A,
  199.      * P10A, P10B, P10C, P5A, P5B, P6A, P6B, P7, P8A, P8B, P8C, P8D,
  200.      * P8E, P8F, P8G, P8H, P9, R, RAND, R0, R4, S, SA, SAX, SAXX, SP,
  201.      * S2, T, T1, T2, T3, T4, T5, T6, T7, X, XC, XE, XEP, XG, XL, XR,
  202.      * X1, X2, X3, X8, Y, YY, YY3, Y2, Y3, BUFFER, SQRTBF, V, SING3,
  203.      * VV, F26
  204.       DIMENSION FDERV(5)
  205.       COMMON /INPUTZ/ A, B, ACCUR, NORM, CHARF, XBREAK, BLEFT, BRIGHT,
  206.      * DBREAK, DEGREE, SMOOTH, LEVEL, EDIST, NBREAK, KNTDIM, NPARDM
  207.       COMMON /RESULZ/ ERROR, KNOTS
  208.       COMMON /KONTRL/ DSCTOL, ERRORI, XLEFT, XRIGHT, BREAK, BOTH,
  209.      * FACTOR, FMESGE, IBREAK, INTERP, LEFT, MAXAUX, MAXKNT, MAXPAR,
  210.      * MAXSTK, NPAR, NSTACK, RIGHT, DISCRD, FATAL, FINISH
  211.       COMMON /COMDIF/ DDTEMP, FINTRP, FLEFT, FRIGHT, XDD, XINTRP,
  212.      * LEFTX, NINTRP, RIGHTX
  213. C
  214.       COMMON /FDATA/ JFUNK, KOUNT, MAXK
  215.       COMMON /FPARS/ PAR2, PAR3A, PAR3B, PAR4, PAR4A, P5A, P5B, P6A,
  216.      * P6B, P7, P8A, P8B, P8C, P8D, P8E, P8F, P8G, P8H, P9, P10A, P10B,
  217.      * P10C
  218. C
  219.       DATA SING3 /.30770707707077/
  220.       DATA BUFFER /1.E-12/
  221.       DATA SQRTBF /1.E-7/
  222. C        DEFINITION OF FUNCTION AT 2700
  223.       F26(VV) = ALOG(2.+VV*VV/(1.+VV))
  224.       F27(V) = SIN(V**2-3.1*EXP(V/(1.+V**2)))*(F26(V)+V)
  225. C
  226. C         FACILITY TO TERMINATE ADAPT TESTING BECAUSE OF EXCESSIVE
  227. C         NUMBER OF FUNCTION EVALUATIONS.
  228.       KOUNT = KOUNT + 1
  229.       IF (KOUNT.NE.MAXK) GO TO 10
  230.       WRITE (6,99999) FMESGE, MAXK
  231.       IF (.NOT.FINISH) FATAL = .TRUE.
  232.       FINISH = .TRUE.
  233.    10 CONTINUE
  234. C
  235. C         PROTECT ARGUMENT X, INITIALIZE DERVS 3,4,5 = 0
  236.       Y = X
  237.       FDERV(3) = 0.0
  238.       FDERV(4) = 0.0
  239.       FDERV(5) = 0.0
  240. C
  241. C         SELECT FUNCTION NUMBER JFUNK
  242.       GO TO (20, 40, 80, 140, 150, 160, 170, 180, 200, 210, 260, 270,
  243.      * 280, 290, 300, 310, 320, 330, 360, 380), JFUNK
  244. C
  245. C        SIMPLE ALBEGRAIC SINGULARITY AT X= 0.0
  246.    20 F = ABS(Y)**.4
  247.       IF (ABS(Y).LE.BUFFER) Y = BUFFER
  248.       FDERV(1) = .4*ABS(Y)**(-.6)*SIGN(1.,Y)
  249.       DO 30 K=2,5
  250.         FK = K
  251.         FDERV(K) = (1.4-FK)*FDERV(K-1)/Y
  252.    30 CONTINUE
  253.       RETURN
  254. C
  255. C         INTERIOR SINGULARITY AT .307707077070770707707077...
  256. C         DIFFERENT EXPONENTS (.481,.63) ON EACH SIDE
  257.    40 T = Y - SING3
  258.       IF (T.LT.0.0) GO TO 60
  259.       F = T**.61
  260.       IF (T.LE.BUFFER) T = BUFFER
  261.       FDERV(1) = .61*F/T
  262.       DO 50 K=2,5
  263.         FDERV(K) = (1.61-FLOAT(K))*FDERV(K-1)/T
  264.    50 CONTINUE
  265.       RETURN
  266.    60 F = (-T)**.481
  267.       FDERV(1) = .481*F/T
  268.       DO 70 K=2,5
  269.         FK = K
  270.         FDERV(K) = (1.481-FK)*FDERV(K-1)/T
  271.    70 CONTINUE
  272.       RETURN
  273. C
  274. C             BROKEN LINE - BREAKS AT 1,2,3.5,4
  275.    80 CONTINUE
  276.       FDERV(2) = 0.0
  277.       IF (Y.LE.1.) GO TO 90
  278.       IF (Y.GE.5.) GO TO 130
  279.       IY = Y + 1.
  280.       GO TO (90, 100, 110, 120, 130), IY
  281.    90 F = 6.*Y + 1.
  282.       FDERV(1) = 6.
  283.       RETURN
  284.   100 F = 7. - (Y-1.)*.5
  285.       FDERV(1) = -.5
  286.       RETURN
  287.   110 F = 6.5 - 1.5*(Y-2.)
  288.       FDERV(1) = -1.5
  289.       RETURN
  290.   120 IF (Y.LT.3.5) GO TO 110
  291.       F = 4.25 - 2.5*(Y-3.5)
  292.       FDERV(1) = -2.5
  293.       RETURN
  294.   130 F = 3. - 3.*(Y-4.)
  295.       FDERV(1) = -3.0
  296.       RETURN
  297. C
  298. C         HUMP AT 2.5 - ONLY 1ST + 2ND DERIVATIVES GIVEN
  299.   140 YY = Y - 2.5
  300.       F = 1./(1.+YY**4)
  301.       YY3 = -4.*YY**3
  302.       FDERV(1) = YY3*F**2
  303.       FDERV(2) = (2.*YY3**2)*F**3 - 12.*YY**2*F**2
  304.       RETURN
  305. C
  306. C           SIMPLE CUBIC
  307.   150 YY = 4.*Y - 2.
  308.       F = YY**3 - YY
  309.       FDERV(1) = 12.*YY**2 - 4.
  310.       FDERV(2) = 96.*YY
  311.       FDERV(3) = 382.
  312.       RETURN
  313. C
  314. C           VARIABLE OSCILLATION - ONLY CORRECT DERIVATIVES ARE 1,2,3
  315.   160 YY = Y**2 - 1.
  316.       F = SIN(YY)
  317.       CS = COS(YY)
  318.       FDERV(1) = 2.*Y*CS
  319.       FDERV(2) = 2.*CS - 4.*(YY+1.)*F
  320.       FDERV(3) = CS*(YY+1.)**3 - 12.*Y*F
  321.       RETURN
  322. C
  323. C           EXPONENTIAL
  324.   170 F = EXP(Y)
  325.       FDERV(1) = F
  326.       FDERV(2) = F
  327.       FDERV(3) = F
  328.       FDERV(4) = F
  329.       FDERV(5) = F
  330.       RETURN
  331. C
  332. C         FIFTH DEGREE POLYNOMIAL PLUS HIGH FREQUENCY SAWTOOTH
  333. C                   ONLY 2 DERIVATIVES GIVEN
  334.   180 F = Y**2*(Y**3-2.) + 6.37
  335.       FD1 = 5.*Y**4 - 4.*Y
  336.       FD2 = 20.*Y**3 - 4.
  337. C         EPS IS FRACTIONAL PART OF 1000*Y
  338.       L = 1000.*Y
  339.       EPS = Y - .001*FLOAT(L)
  340.       L = L - 2*(L/2)
  341.       IF (L.EQ.1) GO TO 190
  342. C             INCREASING PART OF THE SAWTOOTH
  343.       F = F + 500.*EPS**2
  344.       FDERV(1) = FD1 + 1000.*EPS
  345.       FDERV(2) = FD2 + 1000.
  346.       RETURN
  347. C             DECREASING PART OF THE SAWTOOTH
  348.   190 CONTINUE
  349.       F = F + 500.*(1.E-6-EPS**2)
  350.       FDERV(1) = FD1 + 1. - 1000.*EPS
  351.       FDERV(2) = FD2 - 1000.
  352.       RETURN
  353. C
  354. C           SEVENTH DEGREE POLYNOMIAL
  355.   200 F = Y**7 - 2.*Y**6 + 3.1*Y**4 - 6.2*Y
  356.       FDERV(1) = 7.*Y**6 - 12.*Y**5 + 12.4*Y**3 - 6.2
  357.       FDERV(2) = 42.*Y**5 - 60.*Y**4 + 37.2*Y**2
  358.       FDERV(3) = 210.*Y**4 - 240.*Y**3 + 74.4*Y
  359.       FDERV(4) = 840.*Y**3 - 720.*Y**2 + 74.4
  360.       FDERV(5) = 2520.*Y**2 - 1440.*Y
  361.       RETURN
  362. C
  363. C         WILD VARIATION BETWEEN .47 AND .48 - ONLY 3 DERIVATIVES HERE
  364. C         THERE ARE JUMPS AND INFINITIES IN THE SLOPE.
  365.   210 FDERV(2) = 0.
  366.       IF (ABS(Y-.35).GT..13) GO TO 230
  367.       IF (ABS(Y-.45).GT..03) GO TO 240
  368.       IF (ABS(Y-.475).GT..005) GO TO 250
  369. C                SMALL PART - PROB NOT EXACTLY CONTINOUS
  370.       Y2 = -(Y**2-.95*Y+.2256)
  371.       Y3 = ABS(Y2)**.3
  372.       F = 547.5*Y - 259.39 + 26.7*Y3
  373.       IF (Y2.LE.BUFFER) GO TO 220
  374.       FDERV(1) = 547.5 + 8.01*Y3/Y2*(2.*Y-.95)
  375.       FDERV(2) = (-11.214*(2.*Y-.95)**2/Y2+16.02)*Y3/Y2
  376.       FDERV(3) = ((19.0638*(2.*Y-.95)/Y2-44.856)*(2.*Y-.95)-11.214)*Y3/
  377.      * Y2**2
  378.       RETURN
  379. C             SET DERIVATIVES = 0 AT BREAKS
  380.   220 FDERV(1) = 0.0
  381.       FDERV(2) = 0.0
  382.       FDERV(3) = 0.0
  383.       RETURN
  384. C                 MAIN LINEAR PIECE
  385.   230 F = 6.*Y + .55
  386.       FDERV(1) = 6.
  387.       RETURN
  388. C                 NEXT LINEAR PART
  389.   240 F = 25.2*Y - 3.674
  390.       FDERV(1) = 25.2
  391.       RETURN
  392. C                 SMALL LINEAR PART
  393.   250 F = -179.05*Y + 82.111
  394.       FDERV(1) = -179.05
  395.       RETURN
  396. C
  397. C          LOTS OF OSCILLATIONS - FIRST EXAMPLE - ONLY 2 DERIVATIVES
  398.   260 CONTINUE
  399.       X1 = .5*Y**2 + Y
  400.       X2 = .1*Y**3 - 1.
  401.       S = SIN(X1)
  402.       C1 = COS(X1)
  403.       C = COS(X2)
  404.       S2 = SIN(X2)
  405.       R0 = 1./(1.+Y**6)
  406.       R4 = 1./(1.+(Y-4.)**6)
  407.       SP = C1*(Y+1.)
  408.       CP = -S2*(.3*Y**2)
  409.       F = (S+C)*(R0+R4)
  410.       FDERV(1) = -6.*(S+C)*(R0**2*Y**5+R4**2*(Y-4.)**5) +
  411.      * (R0+R4)*(SP+CP)
  412.       FDERV(2) = (S+C)*(-30.*(Y**4*R0**2+(Y-4.)**4*R4**2)+72.*(Y**10*R0*
  413.      * *3+(Y-4.)**10*R4**3)) - 12.*(R0**2*Y**5+(Y-4.)**5*R4**2)*(SP+CP)
  414.      * + (R0+R4)*(C1-(Y+1.)**2*S-.09*Y**4*C-.6*Y*S2)
  415.       RETURN
  416. C
  417. C          LOTS OF OSCILLATIONS - SECOND EXAMPLE - ONLY 3 DERIVATIVES
  418.   270 CONTINUE
  419.       EX = EXP(Y)
  420.       R = 1./(Y+PAR2)
  421.       S = SIN(R)
  422.       C = COS(R)
  423.       DS = -C*R**2
  424.       D2S = -S*R**4 + 2.*C*R**3
  425.       D3S = C*R**6 + 6.*S*R**5 - 6.*C*R**4
  426.       F = EX*S
  427.       FDERV(1) = EX*(S+DS)
  428.       FDERV(2) = EX*(S+2.*DS+D2S)
  429.       FDERV(3) = EX*(S+3.*DS+3.*D2S+D3S)
  430.       RETURN
  431. C
  432. C         END POINT SINGURITIES OF PAR3A,PAR3B - ONLY 3 DERIVATIVES
  433.   280 CONTINUE
  434.       IF (Y.LE.BUFFER) Y = BUFFER
  435.       IF (Y.GE.1.-BUFFER) Y = 1. - BUFFER
  436.       FL = Y**PAR3A
  437.       FR = (1.-Y)**PAR3B
  438.       DL = PAR3A*FL/Y
  439.       DR = -PAR3B*FR/(1.-Y)
  440.       D2L = (PAR3A-1.)*DL/Y
  441.       D2R = -(PAR3B-1.)*DR/(1.-Y)
  442.       D3L = (PAR3A-2.)*D2L/Y
  443.       D3R = -(PAR3B-2.)*D2R/(1.-Y)
  444.       F = FL*FR
  445.       FDERV(1) = FL*DR + DL*FR
  446.       FDERV(2) = FL*D2R + 2.*DL*DR + D2L*FR
  447.       FDERV(3) = FL*D3R + 3.*DL*D2R + 3.*D2L*DR + D3L*FR
  448.       RETURN
  449. C
  450. C         RANDOM NUMBERS ADDED TO SIN(X) - ONLY TWO DERIVATIVES
  451.   290 CONTINUE
  452.       S = SIN(Y)
  453.       C = COS(Y)
  454. C                   QUICKIE PSEUDO-RANDOM NUMBER GENERATOR USED HERE
  455.       RAND = AMOD(4829.*S+C,7.23)/7.23
  456.       F = S + PAR4*(RAND-.5)
  457.       RAND = AMOD(2993.*F+S,3.19)/3.19
  458.       FDERV(1) = C + PAR4*(RAND-.5)
  459.       RAND = AMOD(4901.*RAND+C,8.13)/8.13
  460.       FDERV(2) = -S + PAR4*(RAND-.5)
  461.       RAND = AMOD(6987.*RAND+F*C,4.27)/4.27
  462.       FDERV(3) = -C + PAR4*(RAND-.5)
  463.       RETURN
  464. C
  465. C         LARGE FUNCTION
  466.   300 CONTINUE
  467.       EX = EXP(Y)*P5A
  468.       X8 = Y**8*P5B
  469.       F = EX + X8
  470.       FDERV(1) = EX + 8.*X8/Y
  471.       FDERV(2) = EX + 56.*X8/Y**2
  472.       FDERV(3) = EX + 336.*X8/Y**3
  473.       FDERV(4) = EX + 1680.*Y**4*P5B
  474.       FDERV(5) = EX + 6720.*Y**3*P5B
  475.       RETURN
  476. C
  477. C          SHIFTED FUNCTION   - ONLY 2 DERIVATIVES
  478.   310 CONTINUE
  479.       EX = EXP(Y-P6A)
  480.       X3 = (Y-P6B)**3
  481.       S = SIN(Y)
  482.       C = COS(Y)
  483.       R = 1./(1.+.5*S)
  484.       F = EX + X3*R
  485.       DR = -.5*C*R**2
  486.       D2R = .5*S*R**2 + .5*C**2*R**3
  487.       FDERV(1) = EX + 3.*(Y-P6B)**2*R + X3*DR
  488.       FDERV(2) = EX + 6.*(Y-P6B)*R + 6.*(Y-P6B)**2*DR + X3*D2R
  489.       RETURN
  490. C
  491. C         NUMERICAL DIFFERENTIATON OF ARITHMETIC STATEMENT FUNCTION
  492. C             ONLY  FIRST, 2ND AND THIRD DERIVATIVES
  493.   320 CONTINUE
  494.       F = F27(Y)
  495.       FR = F27(Y+P7)
  496.       FL = F27(Y-P7)
  497.       FRR = F27(Y+2.*P7)
  498.       FLL = F27(Y-2.*P7)
  499.       FDERV(1) = (FR-FL)/P7*.5
  500.       FDERV(2) = (FR-2.*F+FL)/P7**2
  501.       FDERV(3) = (-FLL+2.*FL-2.*FR+FRR)/P7**3*.5
  502.       RETURN
  503. C
  504. C         RATHER COMPLEX FUNCTION - SUM OF 7 DIFFERENT THINGS
  505. C             ONLY TWO DERIVATIVES
  506. C             INFINITE SLOPE AT P8B AND CUSP AT P8E POSSIBLE
  507.   330 CONTINUE
  508. C                   TERM 1
  509.       SA = SQRT(ABS(Y-P8A))
  510.       T1 = EXP(-SA)
  511.       IF (SA.LE.SQRTBF) SA = SQRTBF
  512.       DT1 = -T1*SIGN(1.,Y-P8A)*.5/SA
  513.       DDT1 = T1*(1.+1./SA)*.25/(SA*SA)
  514. C                   TERM 2
  515.       SA = 4.*EXP(-ABS(Y-P8B))
  516.       SAX = Y*SA
  517.       SAXX = Y*SAX
  518.       T2 = 2. - SAXX
  519.       DT2 = -2.*SAX + SAXX*SIGN(1.,Y-P8B)
  520.       DDT2 = -2.*SA + 4.*SAX*SIGN(1.,Y-P8B) - SAXX
  521. C                   TERM 3
  522.       XC = Y - P8C
  523.       T3 = 1./(1.+150.*XC**6)
  524.       DT3 = -900.*XC**5*T3**2
  525.       DDT3 = 1620000.*XC**10*T3**3 - 4500.*XC**4*T3**2
  526. C                   TERM 4
  527.       T4 = 6./(1.+P8D*XC**4)
  528.       DT4 = -2./3.*P8D*XC**3*T4**2
  529.       DDT4 = 8./9.*P8D**2*XC**6*T4**3 - 2.*P8D*XC**2*T4**2
  530. C
  531. C                   TERM 5
  532.       XE = Y - P8E
  533.       XEP = ABS(XE)**P8F
  534.       R = 10./(1.+4.*XE**4)
  535.       T5 = R*XEP
  536.       DR = -1.6*XE**3*R**2
  537.       DDR = 5.12*XE**6*R**3 - 4.8*XE**2*R**2
  538.       DXEP = 0.0
  539.       DDXEP = 0.0
  540.       IF (ABS(XE).LE.1.E-12) GO TO 340
  541.       DXEP = P8F*XEP/XE
  542.       DDXEP = (P8F-1.)*DXEP/XE
  543.   340 DT5 = DR*XEP + R*DXEP
  544.       DDT5 = DDR*XEP + 2.*DR*DXEP + R*DDXEP
  545. C
  546. C                   TERM 6
  547.       XG = ABS(Y-P8G)**2.5
  548.       T6 = -2.*EXP(-XG)
  549.       DT6 = -2.5*T6*SIGN(1.,Y-P8G)*XG**.6
  550.       DDT6 = T6*6.25*(XG**1.2-XG**.2*.6)
  551. C
  552. C                   TERM 7
  553.       S = SIN(4.*Y)
  554.       C = COS(4.*Y)
  555.       R = 1./(.5+ABS(Y-P8H)*5.)
  556.       EX = EXP(Y+1.-P8H)
  557. C             ACTUALLY USING EXP(AMAX1(X-15,0)) - SPLIT CALCULATION
  558.       T7 = S*R
  559.       DT7 = 4.*C*R - S*R**2*SIGN(1.,Y-P8H)*5.
  560.       DDT7 = -16.*S*R - 40.*C*R**2*SIGN(1.,Y-P8H) + 50.*S*R**3
  561.       IF (Y.LT.P8H-1.) GO TO 350
  562. C                        INCLUDE EX TERM
  563.       DDT7 = EX*(DDT7+2.*DT7+T7)
  564.       DT7 = EX*(DT7+T7)
  565.       T7 = EX*T7
  566.   350 CONTINUE
  567. C
  568.       F = T1 + T2 + T3 + T4 + T5 + T6 + T7
  569.       FDERV(1) = DT1 + DT2 + DT3 + DT4 + DT5 + DT6 + DT7
  570.       FDERV(2) = DDT1 + DDT2 + DDT3 + DDT4 + DDT5 + DDT6 + DDT7
  571.       RETURN
  572. C
  573. C         JUMP DISCONTINUITY (EXCEPT FOR P9 = 0.) - ONLY 2 DERIVATIVES
  574.   360 CONTINUE
  575.       IF (Y.GE.P9) GO TO 370
  576. C                     EXPONENTIAL
  577.       F = EXP(Y)
  578.       FDERV(1) = F
  579.       FDERV(2) = F
  580.       FDERV(3) = F
  581.       RETURN
  582.   370 CONTINUE
  583. C                     RATIONAL
  584.       F = 1./(1.+Y**2)
  585.       FDERV(1) = -2.*Y*F**2
  586.       FDERV(2) = 2.*F**2*(4.*Y*F-1.)
  587.       RETURN
  588. C
  589. C                    INFINITE SINGULARITY
  590.   380 CONTINUE
  591.       IF (Y.GT.P10A+BUFFER) GO TO 390
  592.       IF (Y.GT.P10A-BUFFER) GO TO 400
  593. C              LEFT SIDE OF SINGULARITY
  594.       XL = P10A - Y
  595.       F = XL**P10B
  596.       FDERV(1) = -P10B*F/XL
  597.       FDERV(2) = -(P10B-1.)*FDERV(1)/XL
  598.       FDERV(3) = -(P10B-2.)*FDERV(2)/XL
  599.       RETURN
  600.   390 CONTINUE
  601. C              RIGHT SIDE OF SINGULARITY
  602.       XR = Y - P10A
  603.       F = XR**P10C
  604.       FDERV(1) = P10C*F/XR
  605.       FDERV(2) = (P10C-1.)*FDERV(1)/XR
  606.       FDERV(3) = (P10C-2.)*FDERV(2)/XR
  607.       RETURN
  608.   400 CONTINUE
  609. C              AT THE SINGULARITY
  610.       F = 100./(BUFFER**AMAX1(P10B,P10C))
  611.       FDERV(1) = 0.0
  612.       FDERV(2) = F**3
  613.       RETURN
  614. 99999 FORMAT (//2X, 6A4, 14HEXCEEDED LIMIT, I5, 14H ON F(X) EVALS)
  615.       END
  616.       REAL FUNCTION QPOLY(T)
  617. C
  618. C **  THIS IS FUNCTION EXACTLY LIKE PPOLY EXCEPT IT HAS JUST 1 ARGUMENT
  619. C         SO IT CAN BE USED WITH ERRINT BY SUMRY2.
  620. C
  621.       REAL A, ACCUR, B, BLEFT, BRIGHT, CHARF, DDTEMP, DSCTOL, ERROR,
  622.      * ERRORI, FACTOR, FINTRP, FLEFT, FRIGHT, NORM, XBREAK, XDD,
  623.      * XINTRP, XLEFT, XRIGHT
  624.       DIMENSION XBREAK(20), DBREAK(20), BLEFT(20), BRIGHT(20)
  625.       DIMENSION XLEFT(50), XRIGHT(50), FACTOR(12), FMESGE(6)
  626.       DIMENSION DDTEMP(12,12), FINTRP(10), FLEFT(6), FRIGHT(6),
  627.      * XDD(12), XINTRP(10)
  628.       INTEGER BOTH, BREAK, DBREAK, DEGREE, EDIST, FMESGE, RIGHT,
  629.      * RIGHTX, SMOOTH
  630.       LOGICAL DISCRD, FATAL, FINISH
  631.       REAL T, X, ZCOEF, ZKNOTS
  632.       DIMENSION ZKNOTS(100), ZCOEF(100,7)
  633.       COMMON /INPUTZ/ A, B, ACCUR, NORM, CHARF, XBREAK, BLEFT, BRIGHT,
  634.      * DBREAK, DEGREE, SMOOTH, LEVEL, EDIST, NBREAK, KNTDIM, NPARDM
  635.       COMMON /RESULZ/ ERROR, KNOTS
  636.       COMMON /KONTRL/ DSCTOL, ERRORI, XLEFT, XRIGHT, BREAK, BOTH,
  637.      * FACTOR, FMESGE, IBREAK, INTERP, LEFT, MAXAUX, MAXKNT, MAXPAR,
  638.      * MAXSTK, NPAR, NSTACK, RIGHT, DISCRD, FATAL, FINISH
  639.       COMMON /COMDIF/ DDTEMP, FINTRP, FLEFT, FRIGHT, XDD, XINTRP,
  640.      * LEFTX, NINTRP, RIGHTX
  641. C
  642.       COMMON /PARAMS/ ZKNOTS, ZCOEF
  643.       COMMON /SAVEIT/ IKNOT
  644.    10 IF (T.LT.ZKNOTS(IKNOT+1)) GO TO 20
  645.       IKNOT = IKNOT + 1
  646.       IF (IKNOT.LT.KNOTS) GO TO 10
  647.       IKNOT = KNOTS - 1
  648.       GO TO 30
  649.    20 IF (T.GE.ZKNOTS(IKNOT)) GO TO 30
  650.       IKNOT = IKNOT - 1
  651.       IF (IKNOT.LT.1) GO TO 20
  652.       IKNOT = 1
  653.    30 X = T - ZKNOTS(IKNOT)
  654.       QPOLY = ZCOEF(IKNOT,NPAR)
  655.       DO 40 K=1,DEGREE
  656.         KK = NPAR - K
  657.         QPOLY = ZCOEF(IKNOT,KK) + X*QPOLY
  658.    40 CONTINUE
  659.       RETURN
  660.       END
  661.       SUBROUTINE SUMRY2(XKNOTS, COEFS, KDIMEN, NDIMEN)
  662. C
  663. C        =============================================================
  664. C
  665. C **  THIS PROGRAM SUMMARIZES THE PERFORMANCE OF ADAPT
  666. C     IT COMPUTES EXECUTION TIME, INDEPENDENTLY CHECKS THE ACCURACY ,
  667. C         PLOTS F(X) PLUS THE ERROR FOR 100 POINTS.
  668. C     IT USES THE SYSTEM DEPENDENT ROUTINES SECOND AND GRAPH
  669. C         AND THE SUBPROGRAM ERRINT FROM ADAPT.
  670. C
  671.       REAL A, ACCUR, B, BLEFT, BRIGHT, CHARF, DDTEMP, DSCTOL, ERROR,
  672.      * ERRORI, FACTOR, FINTRP, FLEFT, FRIGHT, NORM, XBREAK, XDD,
  673.      * XINTRP, XLEFT, XRIGHT
  674.       DIMENSION XBREAK(20), DBREAK(20), BLEFT(20), BRIGHT(20)
  675.       DIMENSION XLEFT(50), XRIGHT(50), FACTOR(12), FMESGE(6)
  676.       DIMENSION DDTEMP(12,12), FINTRP(10), FLEFT(6), FRIGHT(6),
  677.      * XDD(12), XINTRP(10)
  678.       INTEGER BOTH, BREAK, DBREAK, DEGREE, EDIST, FMESGE, RIGHT,
  679.      * RIGHTX, SMOOTH
  680.       LOGICAL DISCRD, FATAL, FINISH
  681.       REAL XKNOTS(KDIMEN), COEFS(KDIMEN,NDIMEN)
  682.       REAL ABSC(4), DX, DXG, ECHECK, ERRINT, ERRP, P, PPOLY, QPOLY,
  683.      * WGTS(4), XL, XR, ZCOEF, ZKNOTS, F
  684.       REAL XVALU(100), GRAF1(100), GRAF2(100), TIME, TSTART, TSTOP,
  685.      * ZKNOTS(100), ZCOEF(100,7), DUMMY(6)
  686.       EXTERNAL F, QPOLY
  687.       COMMON /INPUTZ/ A, B, ACCUR, NORM, CHARF, XBREAK, BLEFT, BRIGHT,
  688.      * DBREAK, DEGREE, SMOOTH, LEVEL, EDIST, NBREAK, KNTDIM, NPARDM
  689.       COMMON /RESULZ/ ERROR, KNOTS
  690.       COMMON /KONTRL/ DSCTOL, ERRORI, XLEFT, XRIGHT, BREAK, BOTH,
  691.      * FACTOR, FMESGE, IBREAK, INTERP, LEFT, MAXAUX, MAXKNT, MAXPAR,
  692.      * MAXSTK, NPAR, NSTACK, RIGHT, DISCRD, FATAL, FINISH
  693.       COMMON /COMDIF/ DDTEMP, FINTRP, FLEFT, FRIGHT, XDD, XINTRP,
  694.      * LEFTX, NINTRP, RIGHTX
  695. C
  696.       COMMON /PARAMS/ ZKNOTS, ZCOEF
  697.       COMMON /FDATA/ JFUNK, KOUNT, MAXK
  698.       COMMON /INFO/ ECHECK, KFUNK, TIME, TSTART
  699.       DATA ABSC(1) /-.86113631159405/
  700.       DATA ABSC(2) /-.33998104358486/
  701.       DATA ABSC(3) /.33998104358486/
  702.       DATA ABSC(4) /.86113631159405/
  703.       DATA WGTS(1) /.34785484513745/
  704.       DATA WGTS(2) /.65214515486255/
  705.       DATA WGTS(3) /.65214515486255/
  706.       DATA WGTS(4) /.34785484513745/
  707. C
  708.       IF (LEVEL.LE.1) RETURN
  709. C
  710. C                  LEVEL 2 OUTPUT
  711. C
  712. C     ***** CALL CLOCK FOR VALUE OF  TSTOP  IN SECONDS *****
  713. C
  714.       TSTOP = 0.0
  715.       TIME = TSTOP - TSTART
  716.       KFUNK = KOUNT
  717. C
  718. C            VERIFY ERROR BY INDEPENDENT CHECK
  719.       ECHECK = 0.
  720.       IF (KNOTS.GT.100 .OR. NPAR.GT.7) GO TO 40
  721. C             SKIP THIS CHECK
  722. C
  723. C      PUT PARAMETERS IN THE COMMON BLOCK PARAMS FOR QPOLY TO USE
  724.       KNOTS1 = KNOTS - 1
  725.       DO 20 J=1,KNOTS1
  726.         ZKNOTS(J) = XKNOTS(J)
  727.         DO 10 K=1,NPAR
  728.           ZCOEF(J,K) = COEFS(J,K)
  729.    10   CONTINUE
  730.    20 CONTINUE
  731. C
  732. C         BREAK (A,B) INTO HALF AGAIN AS MANY PARTS AS KNOTS
  733.       IPART = (3*KNOTS)/2 + 1
  734.       DX = (B-A)/FLOAT(IPART)
  735.       XL = A
  736.       XR = A + DX
  737.       ECHECK = 0.0
  738. C             COMPUTE ERROR ESTIMATE OVER PART
  739.       P = ABS(NORM)
  740.       DO 30 K=1,IPART
  741.         ERRP = ERRINT(F,QPOLY,XL,XR,ABSC,WGTS)
  742.         XL = XR
  743.         XR = XR + DX
  744. C             UPDATE ERROR
  745.         IF (NORM.EQ.3.) ECHECK = AMAX1(ERRP,ECHECK)
  746.         IF (NORM.NE.3.) ECHECK = (ECHECK**P+ERRP)**(1./P)
  747.    30 CONTINUE
  748.    40 WRITE (6,99999) KFUNK, ACCUR, ERROR, ECHECK
  749.       WRITE (6,99998) TIME
  750. C
  751. C            GRAPHICAL OUTPUT
  752.       RETURN
  753. C         ***** IF GRAPH ROUTINE IS AVAILABLE REMOVE THE ABOVE RETURN
  754. C               AND CALL GRAPH ROUTINE BELOW
  755.       WRITE (6,99997)
  756.       DXG = (B-A)/99.
  757.       DO 50 K=1,100
  758.         XVALU(K) = A + FLOAT(K-1)*DXG
  759. C                  DUMMY IS A DUMMY ARRAY BELOW
  760.         GRAF1(K) = F(XVALU(K),DUMMY)
  761.         GRAF2(K) = GRAF1(K) - PPOLY(XVALU(K),XKNOTS,COEFS,KDIMEN,NDIMEN)
  762.    50 CONTINUE
  763. C     ***** INSERT GRAPH CALL HERE FOR 100 POINTS OF (XVALU,GRAF1,GRAF2)
  764.       RETURN
  765. 99999 FORMAT (/10X, 10HADAPT USED, I5, 28H FUNCTION VALUES FOR ERRORS
  766.      * /10X, 11HSPECIFIED =, E15.5, 21H ESTIMATED BY ADAPT =, E15.5,
  767.      * 21H INDEPENDENT CHECK = , E15.5)
  768. 99998 FORMAT (20X, 11HTIME USED =, F8.5, 9H SECONDS )
  769. 99997 FORMAT (41H1       PLOT OF F(X) AND THE ERROR CURVE )
  770.       END
  771.       SUBROUTINE ADAPT(F, XKNOTS, COEFS, KDIMEN, NDIMEN)
  772. C
  773. C     THIS ALGORITHM COMPUTES A PIECEWISE POLYNOMIAL APPROXIMATION
  774. C  OF SPECIFIED SMOOTHNESS, ACCURACY AND DEGREE.  THE INPUT TO THE
  775. C  COMPUTATION IS
  776. C
  777. C   F      - FUNCTION BEING APPROXIMATED. IT MUST PROVIDE VALUES OF
  778. C            DERIVATIVES UP TO THE ORDER OF SMOOTHNESS SPECIFIED FOR
  779. C            THE APPROXIMATION.  THE CALLING SEQUENCE IS F(X,FDERV) AND
  780. C            FDERV CONTAINS THE DERIVATIVES( SEE CONSTRAINT BELOW)
  781. C   A,B    - THE ENDPOINTS OF THE INTERVAL OF APPROXIMATION
  782. C   ACCUR  - THE ACCURACY REQUIRED FOR THE APPROXIMATION
  783. C   SMOOTH - THE SMOOTHNESS REQUIRED FOR THE APPROXIMATION
  784. C              = 0  MEANS CONTINUOUS
  785. C              = 1  MEANS CONTINUOUS SLOPE
  786. C              = 2  MEANS CONTINUOUS SECOND DERIVATIVE, ETC.
  787. C   DEGREE - THE DEGREE OF THE POLYNOMIAL PIECES.
  788. C            MUST HAVE DEGREE GT 2*SMOOTH
  789. C
  790. C     *****  ***** SECONDARY INPUT - ITEMS WITH DEFAULT VALUES POSSIBLE
  791. C   CHARF  - CHARACTERISTIC LENGTH OF THE FUNCTION F(X). PIECES ARE NOT
  792. C            LONGER THAN THIS LENGTH.
  793. C                   DEFAULT=(B-A) IF DEGREE GT 1, ELSE (B-A)/3
  794. C   NORM   - NORM TO MEASURE THE APPROXIMATION ERROR
  795. C              = 1  L1 APPROXIMATION (LEAST DEVIATIONS)
  796. C              = 2  L2 APPROXIMATION (LEAST SQUARES)
  797. C              = 3  TCHEBYCHEFF (MINIMAX) APPROXIMATION
  798. C              =-P  (NEGATIVE VALUE) GENERAL LP APPROXIMATION
  799. C                   DEFAULT= 2
  800. C   NBREAK - NUMBER OF SPECIAL BREAK POINTS IN THE APPROXIMATION.
  801. C            ASSOCIATED INPUT VARIABLES ARE
  802. C              XBREAK(J)  - LOCATION OF BREAK POINTS
  803. C              DBREAK(J)  - DERIVATIVE BROKEN AT XBREAK
  804. C              BLEFT (J)  - VALUE FROM LEFT FOR DBREAK DERIVATIVE
  805. C              BRIGHT(J)  -   -    -   RIGHT -    -        -
  806. C                   DEFAULT = 0
  807. C   LEVEL  - LEVEL OF OUTPUT FROM ADAPT
  808. C              = -1 NO OUTPUT AT ALL EXCEPT FOR FATAL INPUT ERRORS
  809. C              =  0  ERROR CONDITIONS NOTED, FINAL SUMMARY
  810. C              =  1  PRINT THE APPROXIMATION OUT
  811. C              =  2  DETAILS OF THE COMPUTATION
  812. C              =  3 DEBUG OUTPUT,  = 4 LOTS OF DEBUG OUTPUT
  813. C                    DEFAULT = 0
  814. C   EDIST  - SWITCH TO CHANGE FROM PROPORTIONAL ERROR DISTRIBUTION
  815. C            TO FIXED DISTRIBUTION. THIS IS PRIMARILY OF USE IN
  816. C            APPROXIMATION OF FUNCTIONS WITH SINGULARITIES. ONE SHOULD
  817. C            USE NORM = 1. OR SO IN SUCH CASES
  818. C              = 0  PROPORTIONAL DISTRIBUTION
  819. C              = 1  APPROXIMATE FIXED ERROR DISTRIBUTION
  820. C                   ATTEMPTS TO ACHIEVE SPECIFIED ACCURACY VALUE ACCUR
  821. C              = 2  TRUE FIXED ERROR DISTRIBUTION
  822. C
  823. C     THE OUTPUT OF THE COMPUTATION CONSISTS OF 3 PARTS, EACH RETURNED
  824. C     TO THE USER IN A DIFFERENT WAY. THEY ARE
  825. C
  826. C   XKNOTS,COEFS - ARRAYS DEFINING THE PIECEWISE POLYNOMIAL RESULT.
  827. C   COEFS    XKNOTS(K)  = KNOTS OF THE APPROXIMATION ( K = 1 TO KNOTS)
  828. C                         THE LAST ONE IS RIGHT END POINT OF INTERVAL
  829. C            COEFS(K,N) = COEFFICIENT OF (X - XKNOT(K))**(N-1) IN THE
  830. C                         INTERVAL XKNOT(K) TO XKNOT(K+1)
  831. C                           K = 1 TO KNOTS-1  AND  N = 1 TO DEGREE+1
  832. C            THESE ARRAYS ARE PASSED AS ARGUMENTS SO AS TO USE VARIABLE
  833. C            DIMENSIONS.
  834. C               KDIMEN - DIMENSION USED FOR  XKNOTS IN CALLING PROGRAM
  835. C               NDIMEN - COEFS IS DECLARED  COEFS(KDIMEN,NDIMEN) IN THE
  836. C                        CALLING PROGRAM.
  837. C               ***** NOTE ***** SEVERAL SMALL ARRAYS HERE HAVE FIXED
  838. C                        DIMENSIONS THAT LIMIT DEGREE AND THUS NDIMEN
  839. C                        SHOULD NOT EXCEED THIS LIMIT (CURRENTLY = 12)
  840. C
  841. C   PPOLY  - THE PIECEWISE POLYNOMIAL APPROXIMATING FUNCTION.
  842. C            THIS FUNCTION SUBPROGRAM IS AVAILABLE TO THE USER AT THE
  843. C            COMPLETION OF ADAPT.
  844. C
  845. C   RESULZ - A LABELED COMMON BLOCK WITH  ERROR,KNOTS  IN IT
  846. C              KNOTS - NUMBER OF KNOTS OF PPOLY
  847. C              ERROR - ACCURACY OF PPOLY AS ESTIMATED BY ADAPT
  848. C
  849. C  ********** DIMENSION CONSTRAINTS **********
  850. C     MAXKNT - MAX NUMBER OF KNOTS TAKEN FROM USER VIA KDIMEN
  851. C              ARRAYS WITH THIS DIMENSION
  852. C                   COEFS   XKNOTS
  853. C     MAXPAR - MAX NUMBER OF PARAMETERS PER INTERVAL ( = 12 CURRENTLY )
  854. C                   USER PROVIDED NDIMEN MUST HAVE NDIMEN LE MAXPAR
  855. C                   MUST HAVE  DEGREE + 1 LE MAXPAR
  856. C              ARRAYS WITH THIS DIMENSION (OR RELATED VALUES )
  857. C                D       DDTEMP  FDERVL  FDERVR  FDUMB   FACTOR
  858. C                FINTRP  FLEFT   FRIGHT  POWERS  XTEMP   XINTRP  XDD
  859. C              ***** NOTE ***** MAXPAR ALSO AFFECTS ARGUMENT FDERV
  860. C              OF FUNCTION F.  FDERVL, FDERVR ARE ALSO INVOLVED.
  861. C              SHOULD DECLARE FDERV OF SIZE 6 IN F TO BE SAFE.
  862. C     MAXAUX - MAXIMUM NUMBBER OF AUXILIARY INPUT( = 20 NOW ). ARRAYS
  863. C                 XBREAK  DBREAK  BLEFT   BRIGHT
  864. C     MAXSTK - MAX SIZE OF ACTIVE INTERVAL STACK
  865. C              MIN INTERVAL LENGTH IS 2**(-MAXSTK)*(B-A). ARRAYS
  866. C                 XLEFT   XRIGHT
  867. C
  868. C  **********  PORTABILITY CONSIDERATIONS  **********
  869. C
  870. C     THIS PROGRAM IS IN ANSI STANDARD FORTRAN.  IN ADDITION, IT MEETS
  871. C     ALL THE REQUIREMENTS OF THE BELL LABS PORTABLE FORTRAN -PFORT-
  872. C     EXCEPT ONE.  HOLLERITH CHARACTERS ARE PACKED 4/WORD RATHER THAN
  873. C     1/WORD AS SPECIFIED BY PFORT.
  874. C     NEVERTHELESS, THIS PROGRAM IS AFFECTED IN SEVERAL WAYS BY A
  875. C     CHANGE IN MACHINE WORD LENGTH AND CHANGING TO DOUBLE PRECISION.
  876. C     ***** THIS VERSION IS FOR THE MACHINE WITH THE LONGEST SINGLE
  877. C           PRECISION WORD (CDC).  THE LENGTH OF SOME CONSTANTS IN
  878. C           THE SUBPROGRAM COMPUT EXCEEDS THE CAPACITY OF SOME FORTRAN
  879. C           COMPILERS AND CAN PREVENT COMPILATION.
  880. C  INPUT-OUTPUT -- IS OF THREE TYPES.
  881. C         FATAL ERROR MESSAGES - OCCUR IN SETUP,TAKE,PUT AND TERMIN
  882. C                                THEY CANNOT BE SWITCHED OFF
  883. C         USER AND DEBUGGING OUTPUT - CAN BE SWITCHED OFF
  884. C     THESE INVOLVE MANY FORMATS LIKE  E15.8, F12.8, E22.13, ETC.
  885. C     SOME FORTRAN COMPILERS REQUIRE D-FORMAT FOR DOUBLE PRECISION.
  886. C     SOME DO NOT HANDLE E22.13 ON MACHINES OF SHORTER WORD LENGTH.
  887. C
  888. C     SUMARY PRINTS COEFFICIENTS 5 PER LINE, 6 PER LINE IS BETTER
  889. C     FOR SHORTER WORD LENGTHS.  DOUBLE PRECISION ON MANY MACHINES
  890. C     LIMITS ONE TO 4 PER LINE.
  891. C
  892. C  CONSTANTS -- THE GAUSS WEIGHTS AND ABSCISSA IN COMPUTE ARE GIVEN
  893. C               TO 15 DIGITS IN THE COMMENTS
  894. C
  895. C  DOUBLE PRECISION CONVERSION -- REQUIRES FOUR STEPS
  896. C     1. ALL REAL VARIABLES ARE DECLARED DOUBLE PRECISION. THIS IS
  897. C        DONE BY CHANGING REAL TO DOUBLE PRECISION AS ALL REALS ARE
  898. C        EXPLICITLY DECLARED AND ROOM IS LEFT FOR THIS CHANGE.  REAL
  899. C        VARIABLES APPEAR BEFORE INTEGERS IN ALL COMMON BLOCKS.
  900. C
  901. C     2. ADD D TO CONSTANTS( E.G. 1.0 = 1.0D0 ). ADJUST LENGTH OF
  902. C        GAUSS WEIGHTS IN COMPUTE.
  903. C
  904. C     3. CHANGE ABS,AMAX1,FLOAT  AT MANY PLACES
  905. C
  906. C     4. ADJUST THE INTERVAL STACK SIZE = DIMENSIONS OF XLEFT, XRIGHT
  907. C        AND VALUE OF MAXSTK.  ADJUST THE VALUE BUFFER IN PUT TO BE
  908. C        ABOUT .1E-K FOR A MACHINE WITH K+2 DECIMAL DIGITS.
  909. C
  910.       REAL A, ACCUR, B, BLEFT, BRIGHT, CHARF, DDTEMP, DSCTOL, ERROR,
  911.      * ERRORI, FACTOR, FINTRP, FLEFT, FRIGHT, NORM, XBREAK, XDD,
  912.      * XINTRP, XLEFT, XRIGHT
  913.       DIMENSION XBREAK(20), DBREAK(20), BLEFT(20), BRIGHT(20)
  914.       DIMENSION XLEFT(50), XRIGHT(50), FACTOR(12), FMESGE(6)
  915.       DIMENSION DDTEMP(12,12), FINTRP(10), FLEFT(6), FRIGHT(6),
  916.      * XDD(12), XINTRP(10)
  917.       INTEGER BOTH, BREAK, DBREAK, DEGREE, EDIST, FMESGE, RIGHT,
  918.      * RIGHTX, SMOOTH
  919.       LOGICAL DISCRD, FATAL, FINISH
  920.       REAL XKNOTS(KDIMEN), COEFS(KDIMEN,NDIMEN)
  921.       EXTERNAL F
  922.       REAL F
  923. C
  924.       COMMON /INPUTZ/ A, B, ACCUR, NORM, CHARF, XBREAK, BLEFT, BRIGHT,
  925.      * DBREAK, DEGREE, SMOOTH, LEVEL, EDIST, NBREAK, KNTDIM, NPARDM
  926. C                      KNTDIM - KDIMEN, NAME CHANGED TO PUT IN COMMON
  927. C                      NPARDM - NDIMEN, NAME CHANGED TO PUT IN COMMON
  928.       COMMON /RESULZ/ ERROR, KNOTS
  929. C                      KNOTS = FINAL NO. OF KNOTS, INCLUDES B AS ONE.
  930. C                      ERROR = ESTIMATE OF ERROR ACTUALLY ACHIEVED.
  931.       COMMON /KONTRL/ DSCTOL, ERRORI, XLEFT, XRIGHT, BREAK, BOTH,
  932.      * FACTOR, FMESGE, IBREAK, INTERP, LEFT, MAXAUX, MAXKNT, MAXPAR,
  933.      * MAXSTK, NPAR, NSTACK, RIGHT, DISCRD, FATAL, FINISH
  934. C             KONTRL CONTAINS GENERALLY USEFUL VARIABLES
  935. C                      FATAL  - SWITCH FOR DETECTION OF FATAL ERROR
  936. C                               INCLUDING EXCESSIVE INTERVAL SUBDIVISION
  937. C                               WHICH DOES NOT TERMINATE THE COMPUTATION
  938. C                      FINISH - SWITCH TO TERMINATE ALGORITHM
  939. C                      MAXS   - SEE COMMENTS EARLIER
  940. C                      NSTACK - COUNTER FOR INTERVAL STACK, CONSISTS OF
  941. C                               (XLEFT(J),XRIGHT(J))  J = 1 TO NSTACK
  942. C                      ERRORI - ERROR ESTIMATE FOR TOP INTERVAL
  943. C                      DSCTOL - TOLERANCE TO CHECK DISCARDING INTERVALS
  944. C                      DISCRD - SWITCH TO SIGNAL DISCARD OF TOP INTERVAL
  945. C                      FACTOR - ARRAY OF FACTORIALS
  946. C                      NPAR   - NUMBER OF PAREMETERS = DEGREE + 1
  947. C                      FMESGE - STRING =  **** FATAL ERROR *****
  948. C                      INTERP - NUMBER OF INTERIOR INTERPOLATION POINTS
  949. C                               IN THE NORMAL INTERVAL
  950. C                      IBREAK - COUNTER ON BREAK POINTS
  951. C                      BREAK  - SWITCH FOR BREAK POINT IN TOP INTERVAL
  952. C                               0     = NO BREAK PRESENT
  953. C                               LEFT  = BREAK AT XLEFT(NSTACK)
  954. C                               RIGHT = BREAK AT XRIGHT(NSTACK)
  955. C                               BOTH  = BREAK AT BOTH ENDS
  956.       COMMON /COMDIF/ DDTEMP, FINTRP, FLEFT, FRIGHT, XDD, XINTRP,
  957.      * LEFTX, NINTRP, RIGHTX
  958. C             COMDIF CONTAINS VARIABLES USED ONLY BY COMPUT AND FRIENDS.
  959. C                      NINTRP - NUMBER OF INTERIOR INTERPOLATION POINTS
  960. C                               FOR THE CURRENT INTERVAL
  961. C                      XINTRP - INTERIOR INTERPOLATION POINTS
  962. C                      FINTRP - F VALUES AT XINTRP POINTS
  963. C                      LEFTX  - MULTIPLICITY OF INTERPOLATION AT XLEFT
  964. C                               = NO. OF DERIVATIVES MATCHED AT XLEFT
  965. C                      FLEFT  - VALUES OF F AND ITS DERIVATIVES AT XLEFT
  966. C                      RIGHTX - MULTIPLICITY OF INTERPOLATION AT XRIGHT
  967. C                      FRIGHT - VALUES OF F AND DERIVATIVES AT XRIGHT
  968. C                      DDTEMP - THE ARRAY OF DIVIDED DIFFERENCES
  969. C                      XDD    - THE X VALUES FOR DDTEMP WITH PROPER
  970. C                               MULTIPLICITIES OF XLEFT AND XRIGHT
  971.       COMMON /SAVEIT/ IKNOT
  972. C
  973. C------------------------ MAIN CONTROL PROGRAM -------------------------
  974. C
  975. C   ***** NOTE - ARGUMENTS BELOW ARE FOR READABILITY ONLY    *****
  976. C   *****        EXCEPT FOR F AND XKNOTS,COEFS,KDIMEN,NDIMEN *****
  977. C
  978. C                   CHECK INPUT, INITIAL COMPUTATIONS, PRINT PROBLEM
  979. C
  980.       CALL SETUP(XKNOTS, COEFS, KDIMEN, NDIMEN)
  981. C
  982. C         CHECK FOR FATAL ERROR IN PROBLEM SPECIFICATION
  983.       IF (FATAL) RETURN
  984. C
  985. C                   LOOP OVER PROCESSING OF INTERVALS
  986.    10 CALL TAKE(INTERV)
  987. C
  988.       CALL COMPUT(F, APPROX, INTERV)
  989. C
  990. C                   CHECK FOR DISCARDING INTERVALS
  991.       CALL CHECK(FUNCT, CHARCT)
  992. C
  993. C            PUT NEW INTERVALS ON STACK OR DISCARD, UPDATE STATUS
  994.       CALL PUT(INTERV, XKNOTS, COEFS, KDIMEN, NDIMEN)
  995. C
  996.       CALL TERMIN(TEST, AND, PRINT, XKNOTS, KDIMEN)
  997. C
  998.       IF (.NOT.FINISH) GO TO 10
  999. C
  1000.       CALL SUMARY(XKNOTS, COEFS, KDIMEN, NDIMEN)
  1001.       RETURN
  1002.       END
  1003.       SUBROUTINE CHECK(FUNCT, CHAR)
  1004. C
  1005. C        ===============================================================
  1006. C
  1007. C **  THIS PROGRAM CHECKS FOR DISCARDING INTERVAL, APPLIES VARIOUS
  1008. C     TESTS ABOUT DISCARDING INVOLVING EDIST AND CHARF.
  1009. C
  1010.       REAL A, ACCUR, B, BLEFT, BRIGHT, CHARF, DDTEMP, DSCTOL, ERROR,
  1011.      * ERRORI, FACTOR, FINTRP, FLEFT, FRIGHT, NORM, XBREAK, XDD,
  1012.      * XINTRP, XLEFT, XRIGHT
  1013.       DIMENSION XBREAK(20), DBREAK(20), BLEFT(20), BRIGHT(20)
  1014.       DIMENSION XLEFT(50), XRIGHT(50), FACTOR(12), FMESGE(6)
  1015.       DIMENSION DDTEMP(12,12), FINTRP(10), FLEFT(6), FRIGHT(6),
  1016.      * XDD(12), XINTRP(10)
  1017.       INTEGER BOTH, BREAK, DBREAK, DEGREE, EDIST, FMESGE, RIGHT,
  1018.      * RIGHTX, SMOOTH
  1019.       LOGICAL DISCRD, FATAL, FINISH
  1020.       REAL DTEST, DX
  1021.       COMMON /INPUTZ/ A, B, ACCUR, NORM, CHARF, XBREAK, BLEFT, BRIGHT,
  1022.      * DBREAK, DEGREE, SMOOTH, LEVEL, EDIST, NBREAK, KNTDIM, NPARDM
  1023.       COMMON /RESULZ/ ERROR, KNOTS
  1024.       COMMON /KONTRL/ DSCTOL, ERRORI, XLEFT, XRIGHT, BREAK, BOTH,
  1025.      * FACTOR, FMESGE, IBREAK, INTERP, LEFT, MAXAUX, MAXKNT, MAXPAR,
  1026.      * MAXSTK, NPAR, NSTACK, RIGHT, DISCRD, FATAL, FINISH
  1027.       COMMON /COMDIF/ DDTEMP, FINTRP, FLEFT, FRIGHT, XDD, XINTRP,
  1028.      * LEFTX, NINTRP, RIGHTX
  1029. C
  1030. C
  1031.       DISCRD = .FALSE.
  1032.       DX = XRIGHT(NSTACK) - XLEFT(NSTACK)
  1033. C
  1034. C         COMPUTE DTEST FOR IMPLEMENTING VARIOUS TYPES OF ADAPTIVE APPRX
  1035.       IF (EDIST.EQ.0) DTEST = DX*DSCTOL
  1036. C         FOR THE APPROXIMATE FIXED ERROR DISTRIBUTION TYPE WE ESTIMATE
  1037. C         THE FINAL NUMBER OF KNOTS BY( LIMITING IT A LITTLE )
  1038. C             (NSTACK+KNOTS+2)((B-A)/(XRIGHT-A))
  1039.       IF (EDIST.EQ.1) DTEST = DSCTOL/(FLOAT(NSTACK+KNOTS+2)*AMIN1((B-A)/
  1040.      * (XRIGHT(NSTACK)-A),5.))
  1041.       IF (EDIST.EQ.2 .OR. NORM.EQ.3.) DTEST = DSCTOL
  1042. C
  1043. C             CHECK FOR DISCARD OF INTERVAL
  1044.       IF (ERRORI.LE.DTEST) DISCRD = .TRUE.
  1045. C
  1046. C         CHECK CHARACTERISTIC LENGTH OF FUNCTION
  1047.       IF (DX.GE.CHARF) DISCRD = .FALSE.
  1048.       RETURN
  1049.       END
  1050.       SUBROUTINE COMPUT(F, APPROX, INTERV)
  1051. C
  1052. C        ===============================================================
  1053. C
  1054. C **  THIS PROGRAM COMPUTES THE PIECEWISE POLYNOMIAL APPROXIMATION ON
  1055. C     THE CURRENT INTERVAL. IT ALSO ESTIMATES THE ERROR
  1056. C
  1057.       REAL A, ACCUR, B, BLEFT, BRIGHT, CHARF, DDTEMP, DSCTOL, ERROR,
  1058.      * ERRORI, FACTOR, FINTRP, FLEFT, FRIGHT, NORM, XBREAK, XDD,
  1059.      * XINTRP, XLEFT, XRIGHT
  1060.       DIMENSION XBREAK(20), DBREAK(20), BLEFT(20), BRIGHT(20)
  1061.       DIMENSION XLEFT(50), XRIGHT(50), FACTOR(12), FMESGE(6)
  1062.       DIMENSION DDTEMP(12,12), FINTRP(10), FLEFT(6), FRIGHT(6),
  1063.      * XDD(12), XINTRP(10)
  1064.       INTEGER BOTH, BREAK, DBREAK, DEGREE, EDIST, FMESGE, RIGHT,
  1065.      * RIGHTX, SMOOTH
  1066.       LOGICAL DISCRD, FATAL, FINISH
  1067.       REAL ABSC, DX, ERRINT, F, FDERVL, FDERVR, FDUMB, POLYDD, WGTS
  1068.       DIMENSION ABSC(4), WGTS(4), FDERVL(5), FDERVR(5), FDUMB(5)
  1069.       EXTERNAL F, POLYDD
  1070.       COMMON /INPUTZ/ A, B, ACCUR, NORM, CHARF, XBREAK, BLEFT, BRIGHT,
  1071.      * DBREAK, DEGREE, SMOOTH, LEVEL, EDIST, NBREAK, KNTDIM, NPARDM
  1072.       COMMON /RESULZ/ ERROR, KNOTS
  1073.       COMMON /KONTRL/ DSCTOL, ERRORI, XLEFT, XRIGHT, BREAK, BOTH,
  1074.      * FACTOR, FMESGE, IBREAK, INTERP, LEFT, MAXAUX, MAXKNT, MAXPAR,
  1075.      * MAXSTK, NPAR, NSTACK, RIGHT, DISCRD, FATAL, FINISH
  1076.       COMMON /COMDIF/ DDTEMP, FINTRP, FLEFT, FRIGHT, XDD, XINTRP,
  1077.      * LEFTX, NINTRP, RIGHTX
  1078. C
  1079.       EQUIVALENCE (FLEFT(2),FDERVL(1)), (FRIGHT(2),FDERVR(1))
  1080. C     FIFTEEN DIGIT VALUES FOR THESE GAUSS INTEGRATION CONSTANTS ARE
  1081. C  .861136311594053 .339981043584856 .347854845137454 .652145154862546
  1082.       DATA ABSC(1) /-.86113631159405/
  1083.       DATA ABSC(2) /-.33998104358486/
  1084.       DATA ABSC(3) /.33998104358486/
  1085.       DATA ABSC(4) /.86113631159405/
  1086.       DATA WGTS(1) /.34785484513745/
  1087.       DATA WGTS(2) /.65214515486255/
  1088.       DATA WGTS(3) /.65214515486255/
  1089.       DATA WGTS(4) /.34785484513745/
  1090. C
  1091. C             COMPUTE INTERPOLATION INFORMATION
  1092.       NINTRP = DEGREE - 2*SMOOTH - 1
  1093. C
  1094. C          INCREASE NUMBER OF INTERPOLATION POINTS IF BREAK POINTS ARE
  1095. C          SPECIFIED WITH FEWER DERIVATIVES THAN SMOOTH
  1096.       IF (BREAK.EQ.LEFT .OR. BREAK.EQ.RIGHT) NINTRP = NINTRP + SMOOTH -
  1097.      * DBREAK(IBREAK)
  1098.       IF (BREAK.EQ.BOTH) NINTRP = NINTRP + 2*SMOOTH - DBREAK(IBREAK) -
  1099.      * DBREAK(IBREAK+1)
  1100.       IF (NINTRP.EQ.0) GO TO 20
  1101. C             GENERATE EQUAL SPACED INTERPOLATION POINTS
  1102.       DX = (XRIGHT(NSTACK)-XLEFT(NSTACK))/FLOAT(NINTRP+1)
  1103.       DO 10 J=1,NINTRP
  1104.         XINTRP(J) = XLEFT(NSTACK) + FLOAT(J)*DX
  1105.    10 CONTINUE
  1106. C
  1107. C             GET LEFT AND RIGHT F-VALUES, PUT F-VALUE IN FIRST ELEMENT
  1108. C             OF ARRAYS FLEFT AND FRIGHT.  GET DERIVATIVES BACK AS
  1109. C             OTHER ELEMENTS VIA THE SUBARRAYS FDERVL AND FDERVR.
  1110.    20 FLEFT(1) = F(XLEFT(NSTACK),FDERVL)
  1111.       FRIGHT(1) = F(XRIGHT(NSTACK),FDERVR)
  1112.       LEFTX = SMOOTH + 1
  1113.       RIGHTX = LEFTX
  1114. C            GET F-VALUES AT OTHER INTERPOLATION POINTS, IF ANY
  1115.       IF (NINTRP.EQ.0) GO TO 40
  1116.       DO 30 J=1,NINTRP
  1117.         FINTRP(J) = F(XINTRP(J),FDUMB)
  1118.    30 CONTINUE
  1119. C
  1120. C          CHECK FOR BREAK POINTS, MODIFY VALUES IF NECESSARY
  1121.    40 CONTINUE
  1122.       IF (BREAK.NE.LEFT) GO TO 50
  1123.       LEFTX = DBREAK(IBREAK) + 1
  1124.       FLEFT(LEFTX) = BRIGHT(IBREAK)
  1125.    50 IF (BREAK.NE.RIGHT) GO TO 60
  1126.       RIGHTX = DBREAK(IBREAK) + 1
  1127.       FRIGHT(RIGHTX) = BLEFT(IBREAK)
  1128.    60 IF (BREAK.NE.BOTH) GO TO 70
  1129.       LEFTX = DBREAK(IBREAK) + 1
  1130.       RIGHTX = DBREAK(IBREAK+1) + 1
  1131.       FLEFT(LEFTX) = BRIGHT(IBREAK)
  1132.       FRIGHT(RIGHTX) = BLEFT(IBREAK+1)
  1133.    70 CONTINUE
  1134. C
  1135. C           COMPUTE DIVIDED DIFFERENCES, NEWTON FORM OF POLYNOMIAL
  1136.       CALL NEWTON(LEFTX, RIGHTX, NINTRP)
  1137. C
  1138. C         COMPUTE NORM OF ERROR OF THIS APPROMIMATION USING FOUR PTS
  1139. C         ADD 50 PERCENT FUDGE FACTOR
  1140.       ERRORI = ERRINT(F,POLYDD,XLEFT(NSTACK),XRIGHT(NSTACK),ABSC,WGTS)
  1141.       ERRORI = 1.5*ERRORI
  1142.       RETURN
  1143.       END
  1144.       REAL FUNCTION ERRINT(F, FIT, AAA, BBB, POINTS, WEIGHT)
  1145. C
  1146. C        ===============================================================
  1147. C
  1148. C **  THIS FUNCTION DOES A FOUR POINT INTEGRATION RULE FOR THE
  1149. C     ABSOLUTE VALUE OF THE DIFFERENCE OF TWO FUNCTIONS( F AND FIT )
  1150. C                   ABS( F(X) - FIT(X) )**NORM
  1151. C     THE INTEGRATION USES THE POINTS AND WEIGHTS GIVEN AND SCALED
  1152. C     FROM (-1,1) TO (AAA,BBB)
  1153. C
  1154.       REAL A, ACCUR, B, BLEFT, BRIGHT, CHARF, DDTEMP, DSCTOL, ERROR,
  1155.      * ERRORI, FACTOR, FINTRP, FLEFT, FRIGHT, NORM, XBREAK, XDD,
  1156.      * XINTRP, XLEFT, XRIGHT
  1157.       DIMENSION XBREAK(20), DBREAK(20), BLEFT(20), BRIGHT(20)
  1158.       DIMENSION XLEFT(50), XRIGHT(50), FACTOR(12), FMESGE(6)
  1159.       DIMENSION DDTEMP(12,12), FINTRP(10), FLEFT(6), FRIGHT(6),
  1160.      * XDD(12), XINTRP(10)
  1161.       INTEGER BOTH, BREAK, DBREAK, DEGREE, EDIST, FMESGE, RIGHT,
  1162.      * RIGHTX, SMOOTH
  1163.       LOGICAL DISCRD, FATAL, FINISH
  1164.       REAL AAA, ABMID, BA, BBB, F, FDUMB, FIT, P, PJ, POINTS, WEIGHT
  1165.       REAL ER, F1, F2
  1166.       DIMENSION FDUMB(5), POINTS(4), WEIGHT(4)
  1167.       COMMON /INPUTZ/ A, B, ACCUR, NORM, CHARF, XBREAK, BLEFT, BRIGHT,
  1168.      * DBREAK, DEGREE, SMOOTH, LEVEL, EDIST, NBREAK, KNTDIM, NPARDM
  1169.       COMMON /RESULZ/ ERROR, KNOTS
  1170.       COMMON /KONTRL/ DSCTOL, ERRORI, XLEFT, XRIGHT, BREAK, BOTH,
  1171.      * FACTOR, FMESGE, IBREAK, INTERP, LEFT, MAXAUX, MAXKNT, MAXPAR,
  1172.      * MAXSTK, NPAR, NSTACK, RIGHT, DISCRD, FATAL, FINISH
  1173.       COMMON /COMDIF/ DDTEMP, FINTRP, FLEFT, FRIGHT, XDD, XINTRP,
  1174.      * LEFTX, NINTRP, RIGHTX
  1175. C
  1176. C         COMPUTE MIDPOINT = ABMID AND HALF LENGTH = BA OF INTERVAL
  1177.       ABMID = (AAA+BBB)/2.
  1178.       BA = (BBB-AAA)/2.
  1179.       PJ = ABMID + BA*POINTS(1)
  1180. C
  1181. C         TEST FOR TCHEBYCHEFF (MINIMAX) NORM WHICH USES SPECIAL CODE
  1182.       IF (NORM.EQ.3.) GO TO 20
  1183. C
  1184. C         HAVE GENERAL LP NORM OR LEAST SQUARES OR LEAST DEVIATIONS
  1185.       P = ABS(NORM)
  1186. C              INITIALIZE THE QUADRATURE RULE
  1187.       ERRINT = ABS(F(PJ,FDUMB)-FIT(PJ))**P*WEIGHT(1)
  1188. C              LOOP THROUGH REMAINING POINTS
  1189.       DO 10 J=2,4
  1190.         PJ = ABMID + BA*POINTS(J)
  1191. C             DEBUG DEBUG DEBUG START OF THE COMPUTATION
  1192.         F1 = F(PJ,FDUMB)
  1193.         F2 = FIT(PJ)
  1194.         ER = ABS(F1-F2)**P
  1195.         IF (LEVEL.GE.4 .AND. (KNOTS.EQ.1 .OR. FINISH)) WRITE (6,99999)
  1196.      *   PJ, F1, F2, ER
  1197.         ERRINT = ERRINT + ABS(F(PJ,FDUMB)-FIT(PJ))**P*WEIGHT(J)
  1198.    10 CONTINUE
  1199.       ERRINT = ERRINT*BA
  1200.       GO TO 40
  1201. C
  1202. C         TCHEBYCHEFF NORM
  1203.    20 CONTINUE
  1204. C             FIND MAX ERROR ON POINTS
  1205. C               INITIALIZE
  1206.       ERRINT = ABS(F(PJ,FDUMB)-FIT(PJ))
  1207. C              LOOP THROUGH THE REMAINING POINTS
  1208.       DO 30 J=2,4
  1209.         PJ = ABMID + BA*POINTS(J)
  1210.         ERRINT = AMAX1(ERRINT,ABS(F(PJ,FDUMB)-FIT(PJ)))
  1211.    30 CONTINUE
  1212.    40 CONTINUE
  1213. C             DEBUG  DEBUG  DEBUG  DEBUG
  1214.       IF (LEVEL.GE.3) WRITE (6,99998) ERRINT, AAA, BBB
  1215.       RETURN
  1216. 99999 FORMAT (5(3H --), 31HDEBUG ERROR CURVE,PJ,F1,F2,ER =, 4F15.8)
  1217. 99998 FORMAT (15X, 9HERRINT = , F20.15, 4H ON , 2F15.8)
  1218.       END
  1219.       SUBROUTINE NEWTON(NL, NR, NI)
  1220. C
  1221. C        ===============================================================
  1222. C
  1223. C **  THIS PROGRAM COMPUTES THE DIVIDED DIFFERENCES ARRAY AS FOLLOWS
  1224. C         NL COALESCED POINTS ON LEFT   - DERIV VALUES IN FLEFT
  1225. C         NR COALESCED POINTS ON RIGHT  -   -      -    - FRIGHT
  1226. C         NI DISTINCT  POINTS INBETWEEN - FNCTN    -    - FINTRP
  1227. C
  1228. C     THE POINTS ARE ORDERED XL = XLEFT (NSTACK)
  1229. C                            XR = XRIGHT(NSTACK)
  1230. C                            XINTRP ARRAY
  1231. C
  1232. C         LAYOUT OF THE DDTEMP DIVIDED DIFFERENCE ARRAY
  1233. C
  1234. C     NL=6    LLLLLL****II
  1235. C     NR=4    LLLLL****II     L = FIRST TRIANGLE
  1236. C     NI=2    LLLL****II
  1237. C             LLL****II       R = SECOND TRIANGLE
  1238. C             LL****II
  1239. C             L****II         * = FILL BETWEEN TRIANGLES
  1240. C             RRRRII
  1241. C             RRRII           I = COMPLETION FOR INTERPOLATION POINTS
  1242. C             RRII
  1243. C             RII       IDIF = HORIZONTAL COORD. = DIFFERENCE ORDER
  1244. C             II        IPT  = VERTICAL COORD. ASSOCIATED WITH POINTS
  1245. C             I
  1246. C
  1247. C
  1248.       REAL A, ACCUR, B, BLEFT, BRIGHT, CHARF, DDTEMP, DSCTOL, ERROR,
  1249.      * ERRORI, FACTOR, FINTRP, FLEFT, FRIGHT, NORM, XBREAK, XDD,
  1250.      * XINTRP, XLEFT, XRIGHT
  1251.       DIMENSION XBREAK(20), DBREAK(20), BLEFT(20), BRIGHT(20)
  1252.       DIMENSION XLEFT(50), XRIGHT(50), FACTOR(12), FMESGE(6)
  1253.       DIMENSION DDTEMP(12,12), FINTRP(10), FLEFT(6), FRIGHT(6),
  1254.      * XDD(12), XINTRP(10)
  1255.       INTEGER BOTH, BREAK, DBREAK, DEGREE, EDIST, FMESGE, RIGHT,
  1256.      * RIGHTX, SMOOTH
  1257.       LOGICAL DISCRD, FATAL, FINISH
  1258.       REAL DIFFF, DIFFX
  1259.       COMMON /INPUTZ/ A, B, ACCUR, NORM, CHARF, XBREAK, BLEFT, BRIGHT,
  1260.      * DBREAK, DEGREE, SMOOTH, LEVEL, EDIST, NBREAK, KNTDIM, NPARDM
  1261.       COMMON /RESULZ/ ERROR, KNOTS
  1262.       COMMON /KONTRL/ DSCTOL, ERRORI, XLEFT, XRIGHT, BREAK, BOTH,
  1263.      * FACTOR, FMESGE, IBREAK, INTERP, LEFT, MAXAUX, MAXKNT, MAXPAR,
  1264.      * MAXSTK, NPAR, NSTACK, RIGHT, DISCRD, FATAL, FINISH
  1265.       COMMON /COMDIF/ DDTEMP, FINTRP, FLEFT, FRIGHT, XDD, XINTRP,
  1266.      * LEFTX, NINTRP, RIGHTX
  1267. C
  1268. C                 MAIN CALCULATION OF DIVIDED DIFFERENCES
  1269. C         DEFINE A FEW SHORT CONSTANTS
  1270.       NL1 = NL - 1
  1271.       NL2 = NL + 1
  1272.       NR1 = NR - 1
  1273.       NR2 = NR + 1
  1274.       NRL = NR + NL
  1275. C
  1276. C         PUT X-VALUES IN A SINGLE ARRAY WITH NDDX = NL+NR+NI POINTS
  1277.       DO 10 NDDX=1,NL
  1278.         XDD(NDDX) = XLEFT(NSTACK)
  1279.    10 CONTINUE
  1280.       NDDX = NL
  1281.       DO 20 K=1,NR
  1282.         NDDX = NDDX + 1
  1283.         XDD(NDDX) = XRIGHT(NSTACK)
  1284.    20 CONTINUE
  1285. C            CHECK IF THERE ARE ANY INTERPOLATION POINTS TO ADD TO XDD
  1286.       IF (NI.EQ.0) GO TO 40
  1287.       DO 30 K=1,NI
  1288.         NDDX = NDDX + 1
  1289.         XDD(NDDX) = XINTRP(K)
  1290.    30 CONTINUE
  1291. C
  1292. C           FILL BORDER OF FIRST TRIANGLE - SIZE NL.
  1293.    40 CONTINUE
  1294. C         TOP BORDER
  1295.       DO 50 IDIF=1,NL
  1296.         DDTEMP(IDIF,1) = FLEFT(IDIF)/FACTOR(IDIF)
  1297.    50 CONTINUE
  1298.       IF (NL1.EQ.0) GO TO 70
  1299. C                   BOTTOM BORDER
  1300.       DO 60 IDIF=1,NL1
  1301.         IPT = NL2 - IDIF
  1302.         DDTEMP(IDIF,IPT) = DDTEMP(IDIF,1)
  1303.    60 CONTINUE
  1304. C
  1305. C          FILL BORDER OF SECOND TRIANGLE - SIZE NR
  1306.    70 CONTINUE
  1307. C         TOP BORDER
  1308.       DO 80 IDIF=1,NR
  1309.         DDTEMP(IDIF,NL2) = FRIGHT(IDIF)/FACTOR(IDIF)
  1310.    80 CONTINUE
  1311.       IF (NRL.EQ.NL2) GO TO 100
  1312. C                   BOTTOM BORDER
  1313.       DO 90 IDIF=1,NR1
  1314.         IPT = NRL + 1 - IDIF
  1315.         DDTEMP(IDIF,IPT) = DDTEMP(IDIF,NL2)
  1316.    90 CONTINUE
  1317. C
  1318. C           FILL PARALLOGRAM BETWEEN THE TWO TRIANGLES JUST FILLED
  1319. C          FILL ENTRIES PARALLEL TO BOTTOM OF FIRST TRIANGLE
  1320.   100 CONTINUE
  1321. C
  1322. C         LOOP STEPPING ALONG TOP SIDE OF SECOND TRIANGLE
  1323.       DO 120 J=2,NR2
  1324.         IDIF = J
  1325. C             LOOP STEPPING PARALLEL TO BOTTOM SIDE OF FIRST TRIANGLE
  1326.         DO 110 K=2,NL2
  1327.           IPT = NL + 2 - K
  1328.           DIFFF = DDTEMP(IDIF-1,IPT+1) - DDTEMP(IDIF-1,IPT)
  1329.           IPT2 = IPT - 1 + IDIF
  1330.           DIFFX = XDD(IPT2) - XDD(IPT)
  1331.           DDTEMP(IDIF,IPT) = DIFFF/DIFFX
  1332.           IDIF = IDIF + 1
  1333.   110   CONTINUE
  1334.   120 CONTINUE
  1335. C             DEBUG  DEBUG  DEBUG  DEBUG
  1336.       IF (LEVEL.GE.4 .AND. KNOTS.LE.1) WRITE (6,99999) NR2, NL2, IDIF,
  1337.      * IPT, DIFFF, DIFFX
  1338. C
  1339. C         FILL IN BOTTOM DIAGONALS FOR INTERPOLATION POINTS, IF ANY
  1340.       IF (NI.EQ.0) GO TO 150
  1341. C         LOOP THROUGH THE INTERPOLATATION POINTS
  1342.       DO 140 J=1,NI
  1343.         IDIF = 2
  1344.         NRLJ = NRL + J
  1345.         DDTEMP(1,NRLJ) = FINTRP(J)
  1346. C         LOOP THROUGH THE DIFFERENCES (IDIF INDEX)
  1347.         NRLJ1 = NRLJ - 1
  1348.         DO 130 K=1,NRLJ1
  1349.           IPT = NRLJ - K
  1350.           DIFFF = DDTEMP(IDIF-1,IPT+1) - DDTEMP(IDIF-1,IPT)
  1351.           DIFFX = XDD(NRLJ) - XDD(IPT)
  1352.           DDTEMP(IDIF,IPT) = DIFFF/DIFFX
  1353.           IDIF = IDIF + 1
  1354.   130   CONTINUE
  1355.   140 CONTINUE
  1356.   150 CONTINUE
  1357.       RETURN
  1358. 99999 FORMAT (9X, 30HNR2,NL2,IDIF,IPT,DIFFF,DIFFX =, 4I3, 2F12.6)
  1359.       END
  1360.       REAL FUNCTION POLYDD(X)
  1361. C
  1362. C        ===============================================================
  1363. C
  1364. C **  THIS FUNCTION EVALUATES THE CURRENT POLYNOMIAL PIECE REPRESENTED
  1365. C     BY THE DIVIDED DIFFERENCES DDTEMP ON THE POINT SET XDD.
  1366. C
  1367.       REAL A, ACCUR, B, BLEFT, BRIGHT, CHARF, DDTEMP, DSCTOL, ERROR,
  1368.      * ERRORI, FACTOR, FINTRP, FLEFT, FRIGHT, NORM, XBREAK, XDD,
  1369.      * XINTRP, XLEFT, XRIGHT
  1370.       DIMENSION XBREAK(20), DBREAK(20), BLEFT(20), BRIGHT(20)
  1371.       DIMENSION XLEFT(50), XRIGHT(50), FACTOR(12), FMESGE(6)
  1372.       DIMENSION DDTEMP(12,12), FINTRP(10), FLEFT(6), FRIGHT(6),
  1373.      * XDD(12), XINTRP(10)
  1374.       INTEGER BOTH, BREAK, DBREAK, DEGREE, EDIST, FMESGE, RIGHT,
  1375.      * RIGHTX, SMOOTH
  1376.       LOGICAL DISCRD, FATAL, FINISH
  1377.       REAL X
  1378.       COMMON /INPUTZ/ A, B, ACCUR, NORM, CHARF, XBREAK, BLEFT, BRIGHT,
  1379.      * DBREAK, DEGREE, SMOOTH, LEVEL, EDIST, NBREAK, KNTDIM, NPARDM
  1380.       COMMON /RESULZ/ ERROR, KNOTS
  1381.       COMMON /KONTRL/ DSCTOL, ERRORI, XLEFT, XRIGHT, BREAK, BOTH,
  1382.      * FACTOR, FMESGE, IBREAK, INTERP, LEFT, MAXAUX, MAXKNT, MAXPAR,
  1383.      * MAXSTK, NPAR, NSTACK, RIGHT, DISCRD, FATAL, FINISH
  1384.       COMMON /COMDIF/ DDTEMP, FINTRP, FLEFT, FRIGHT, XDD, XINTRP,
  1385.      * LEFTX, NINTRP, RIGHTX
  1386. C
  1387. C
  1388.       POLYDD = DDTEMP(DEGREE+1,1)
  1389.       DO 10 K=1,DEGREE
  1390.         J = DEGREE + 1 - K
  1391.         POLYDD = DDTEMP(J,1) + (X-XDD(J))*POLYDD
  1392.    10 CONTINUE
  1393.       RETURN
  1394.       END
  1395.       SUBROUTINE PPFIT4(F, XLFT, XRGT, EPSLN, NPIECE, ERREST, XKNOTS,
  1396.      * COEFS, KDIMEN, NDIMEN, NDEG, NSMTH, EMEAS, LPRNT, FOSCIL, ATYPE,
  1397.      * KBREAK, BRAKPT, KDERVB, VALLFT, VALRGT)
  1398. C
  1399. C        ===============================================================
  1400. C
  1401. C **  THIS SET OF FOUR CONTROL PROGRAMS SET VARYING NUMBERS OF DEFAULT
  1402. C     VALUES FOR ARGUMENTS.  IT USES ENTRY STATEMENTS WHICH ARE DONE
  1403. C     DIFFERENTLY BY VARIOUS FORTRANS.  FOR THIS REASON ENTRY STATEMENTS
  1404. C     ARE ONLY INDICATED BY COMMENT CARDS.  WRITING FOUR SEPARATE ROU-
  1405. C     TINES APPROXIMATELY TRIPLES THE LENGTH OF THE TOTAL CODE.
  1406. C
  1407. C     THE FOLLOWING TABULATES THE INTERNAL AND EXTERNAL NAMES OF THE
  1408. C     ARGUMENTS ALONG WITH THEIR DEFAULT VALUES FOR THE VARIOUS PPFIT.
  1409. C     NOTE THAT THIS ALLOWS THE ARGUMENTS TO BE PUT INTO A COMMON BLOCK
  1410. C     AND AVOIDS LONG ARGUMENT LISTS WITHIN ADAPT ITSELF.
  1411. C
  1412. C
  1413. C    INTERNAL   EXTERNAL     DEFAULT VALUE    SET BY PPFIT NUMBER
  1414. C     F          F              NONE
  1415. C     A,B        A,B            NONE
  1416. C     ACCUR      EPSLN          NONE
  1417. C     KNOTS      NPIECE         OUTPUT
  1418. C     ERROR      ERREST         OUTPUT
  1419. C     XKNOTS     XKNOTS         OUTPUT
  1420. C     COEFS      COEFS          OUTPUT
  1421. C     KDIMEN     KDIMEN         NONE
  1422. C     NDIMEN     NDIMEN         NONE
  1423. C     DEGREE     NDEG           3                1
  1424. C     SMOOTH     NSMTH          0                1
  1425. C     NORM       EMEAS          3                1 2
  1426. C     LEVEL      LPRNT          1                1 2
  1427. C     CHARF      FOSCIL         VARIABLE         1 2
  1428. C     EDIST      ATYPE          1                1 2 3
  1429. C     NBREAK     KBREAK         0                1 2 3
  1430. C     XBREAK     BRAKPT         -
  1431. C     DBREAK     KDERVB         -
  1432. C     BLEFT      VALLFT         -
  1433. C     BRIGHT     VALRGT         -
  1434. C
  1435.       REAL A, ACCUR, B, BLEFT, BRIGHT, CHARF, DDTEMP, DSCTOL, ERROR,
  1436.      * ERRORI, FACTOR, FINTRP, FLEFT, FRIGHT, NORM, XBREAK, XDD,
  1437.      * XINTRP, XLEFT, XRIGHT
  1438.       DIMENSION XBREAK(20), DBREAK(20), BLEFT(20), BRIGHT(20)
  1439.       DIMENSION XLEFT(50), XRIGHT(50), FACTOR(12), FMESGE(6)
  1440.       DIMENSION DDTEMP(12,12), FINTRP(10), FLEFT(6), FRIGHT(6),
  1441.      * XDD(12), XINTRP(10)
  1442.       INTEGER BOTH, BREAK, DBREAK, DEGREE, EDIST, FMESGE, RIGHT,
  1443.      * RIGHTX, SMOOTH
  1444.       LOGICAL DISCRD, FATAL, FINISH
  1445.       REAL XKNOTS(KDIMEN), COEFS(KDIMEN,NDIMEN)
  1446.       REAL BRAKPT, EMEAS, EPSLN, ERREST, FOSCIL, F, VALLFT, VALRGT,
  1447.      * XLFT, XRGT
  1448.       DIMENSION BRAKPT(20), KDERVB(20), VALLFT(20), VALRGT(20)
  1449.       INTEGER ATYPE
  1450.       EXTERNAL F
  1451.       COMMON /INPUTZ/ A, B, ACCUR, NORM, CHARF, XBREAK, BLEFT, BRIGHT,
  1452.      * DBREAK, DEGREE, SMOOTH, LEVEL, EDIST, NBREAK, KNTDIM, NPARDM
  1453.       COMMON /RESULZ/ ERROR, KNOTS
  1454.       COMMON /KONTRL/ DSCTOL, ERRORI, XLEFT, XRIGHT, BREAK, BOTH,
  1455.      * FACTOR, FMESGE, IBREAK, INTERP, LEFT, MAXAUX, MAXKNT, MAXPAR,
  1456.      * MAXSTK, NPAR, NSTACK, RIGHT, DISCRD, FATAL, FINISH
  1457.       COMMON /COMDIF/ DDTEMP, FINTRP, FLEFT, FRIGHT, XDD, XINTRP,
  1458.      * LEFTX, NINTRP, RIGHTX
  1459.       COMMON /SAVEIT/ IKNOT
  1460. C
  1461. C
  1462. C             SKIP ALL DEFAULT SETTING FOR PPFIT4
  1463.       GO TO 10
  1464. C
  1465. C                 SET DEFAULTS FOR PPFIT1
  1466. C
  1467. C         ****** SINCE ENTRY IS NON-STANDARD  ******
  1468. C         ****** THE NATURAL WAY TO IMPLEMENT ******
  1469. C         ****** THE OTHER PPFITS IS ONLY     ******
  1470. C         ****** INDICATED IN THE COMMENTS    ******
  1471. C
  1472. C     ENTRY PPFIT1
  1473. C     ENTRY PPFIT1(F,XLFT,XRGT,EPSLN,NPIECE,ERREST,XKNOTS,COEFS,
  1474. C    A             KDIMEN,NDIMEN)
  1475.       NDEG = 3
  1476.       NSMTH = 0
  1477. C
  1478. C                 SET DEFAULTS FOR PPFIT2
  1479. C     ENTRY PPFIT2
  1480. C     ENTRY PPFIT2(F,XLFT,XRGT,EPSLN,NPIECE,ERREST,XKNOTS,COEFS,
  1481. C    A             KDIMEN,NDIMEN,NDEG,NSMTH)
  1482.       EMEAS = 3.0
  1483.       LPRNT = 1
  1484.       FOSCIL = B - A
  1485.       IF (NDEG.EQ.2) FOSCIL = .5*FOSCIL
  1486.       IF (NDEG.LE.1) FOSCIL = (B-A)/3.
  1487. C
  1488. C                 SET DEFAULTS FOR PPFIT3
  1489. C     ENTRY PPFIT3
  1490. C     ENTRY PPFIT3(F,XLFT,XRGT,EPSLN,NPIECE,ERREST,XKNOTS,COEFS,
  1491. C    A             KDIMEN,NDIMEN,NDEG,NSMTH,EMEAS,LPRNT,FOSCIL)
  1492.       ATYPE = 1
  1493.       KSING = 0
  1494.       KBREAK = 0
  1495. C
  1496. C         PUT INPUT INTO COMMON INPUTZ BY CHANGING TO INTERNAL NAMES
  1497.    10 A = XLFT
  1498.       B = XRGT
  1499.       ACCUR = EPSLN
  1500.       DEGREE = NDEG
  1501.       SMOOTH = NSMTH
  1502.       NORM = EMEAS
  1503.       LEVEL = LPRNT
  1504.       CHARF = FOSCIL
  1505.       EDIST = ATYPE
  1506.       NBREAK = KBREAK
  1507.       IF (NBREAK.LE.0 .OR. NBREAK.GE.21) GO TO 30
  1508.       DO 20 K=1,NBREAK
  1509.         XBREAK(K) = BRAKPT(K)
  1510.         DBREAK(K) = KDERVB(K)
  1511.         BLEFT(K) = VALLFT(K)
  1512.         BRIGHT(K) = VALRGT(K)
  1513.    20 CONTINUE
  1514.    30 CONTINUE
  1515. C
  1516.       CALL ADAPT(F, XKNOTS, COEFS, KDIMEN, NDIMEN)
  1517.       NPIECE = KNOTS
  1518.       ERREST = ERROR
  1519.       RETURN
  1520.       END
  1521.       REAL FUNCTION PPOLY(T, XKNOTS, COEFS, KDIMEN, NDIMEN)
  1522. C
  1523. C        ===============================================================
  1524. C
  1525. C **  THIS FUNCTION EVALUATES THE PIECEWISE POLYNOMIAL APPROXIMATION
  1526. C     COMPUTED IN ADAPT
  1527. C
  1528.       REAL A, ACCUR, B, BLEFT, BRIGHT, CHARF, DDTEMP, DSCTOL, ERROR,
  1529.      * ERRORI, FACTOR, FINTRP, FLEFT, FRIGHT, NORM, XBREAK, XDD,
  1530.      * XINTRP, XLEFT, XRIGHT
  1531.       DIMENSION XBREAK(20), DBREAK(20), BLEFT(20), BRIGHT(20)
  1532.       DIMENSION XLEFT(50), XRIGHT(50), FACTOR(12), FMESGE(6)
  1533.       DIMENSION DDTEMP(12,12), FINTRP(10), FLEFT(6), FRIGHT(6),
  1534.      * XDD(12), XINTRP(10)
  1535.       INTEGER BOTH, BREAK, DBREAK, DEGREE, EDIST, FMESGE, RIGHT,
  1536.      * RIGHTX, SMOOTH
  1537.       LOGICAL DISCRD, FATAL, FINISH
  1538.       REAL XKNOTS(KDIMEN), COEFS(KDIMEN,NDIMEN)
  1539.       REAL T, X
  1540.       COMMON /INPUTZ/ A, B, ACCUR, NORM, CHARF, XBREAK, BLEFT, BRIGHT,
  1541.      * DBREAK, DEGREE, SMOOTH, LEVEL, EDIST, NBREAK, KNTDIM, NPARDM
  1542.       COMMON /RESULZ/ ERROR, KNOTS
  1543.       COMMON /KONTRL/ DSCTOL, ERRORI, XLEFT, XRIGHT, BREAK, BOTH,
  1544.      * FACTOR, FMESGE, IBREAK, INTERP, LEFT, MAXAUX, MAXKNT, MAXPAR,
  1545.      * MAXSTK, NPAR, NSTACK, RIGHT, DISCRD, FATAL, FINISH
  1546.       COMMON /COMDIF/ DDTEMP, FINTRP, FLEFT, FRIGHT, XDD, XINTRP,
  1547.      * LEFTX, NINTRP, RIGHTX
  1548. C
  1549.       COMMON /SAVEIT/ IKNOT
  1550. C           START SEARCH FOR RIGHT INTERVAL FROM POINT OF LAST
  1551. C           EVALUATION OF PPOLY
  1552. C
  1553. C           FORWARD SEARCHING LOOP
  1554.    10 IF (T.LT.XKNOTS(IKNOT+1)) GO TO 20
  1555.       IKNOT = IKNOT + 1
  1556.       IF (IKNOT.LT.KNOTS) GO TO 10
  1557.       IKNOT = KNOTS - 1
  1558. C
  1559. C                REACHED RIGHT END OF INTERVAL
  1560.       GO TO 30
  1561. C
  1562. C           BACKWARD SEARCHING LOOP
  1563.    20 IF (T.GE.XKNOTS(IKNOT)) GO TO 30
  1564.       IKNOT = IKNOT - 1
  1565.       IF (IKNOT.GT.1) GO TO 20
  1566.       IKNOT = 1
  1567. C
  1568. C                REACHED LEFT END OF INTERVAL
  1569. C
  1570. C                EVALUATE FROM POWERS BASED AT XKNOT(IKNOT)
  1571. C                USE NESTED MULTIPLICATION
  1572.    30 X = T - XKNOTS(IKNOT)
  1573.       PPOLY = COEFS(IKNOT,NPAR)
  1574.       DO 40 K=1,DEGREE
  1575.         KK = NPAR - K
  1576.         PPOLY = COEFS(IKNOT,KK) + X*PPOLY
  1577.    40 CONTINUE
  1578.       RETURN
  1579.       END
  1580.       SUBROUTINE PTRANS(D, POWERS)
  1581. C
  1582. C        ===============================================================
  1583. C
  1584. C **  THIS PROGRAM CONVERTS POLYNOMIAL REPRESENTATION FROM DIVIDED
  1585. C     DIFFERENCE TO POWER FORM.  THERE ARE COALESCED POINTS ON EACH
  1586. C     END OF THE INTERVAL (XL,XR) = (XLEFT(NSTACK),XRIGHT(NSTACK)).
  1587. C     THE NUMBER COALESCED AT EACH END IS LEFTX AND RIGHTX.
  1588. C     AND THERE ARE NINTRP OTHER PTS XINTRP(K)  INBETWEEN THEM.
  1589. C     SEE SUBROUTINE NEWTON FOR MORE DETAILS
  1590. C
  1591.       REAL A, ACCUR, B, BLEFT, BRIGHT, CHARF, DDTEMP, DSCTOL, ERROR,
  1592.      * ERRORI, FACTOR, FINTRP, FLEFT, FRIGHT, NORM, XBREAK, XDD,
  1593.      * XINTRP, XLEFT, XRIGHT
  1594.       DIMENSION XBREAK(20), DBREAK(20), BLEFT(20), BRIGHT(20)
  1595.       DIMENSION XLEFT(50), XRIGHT(50), FACTOR(12), FMESGE(6)
  1596.       DIMENSION DDTEMP(12,12), FINTRP(10), FLEFT(6), FRIGHT(6),
  1597.      * XDD(12), XINTRP(10)
  1598.       INTEGER BOTH, BREAK, DBREAK, DEGREE, EDIST, FMESGE, RIGHT,
  1599.      * RIGHTX, SMOOTH
  1600.       LOGICAL DISCRD, FATAL, FINISH
  1601.       REAL D, POWERS, SHIFT, XL, XR, XTEMP
  1602.       DIMENSION D(12,12), POWERS(12), XTEMP(12)
  1603.       COMMON /INPUTZ/ A, B, ACCUR, NORM, CHARF, XBREAK, BLEFT, BRIGHT,
  1604.      * DBREAK, DEGREE, SMOOTH, LEVEL, EDIST, NBREAK, KNTDIM, NPARDM
  1605.       COMMON /RESULZ/ ERROR, KNOTS
  1606.       COMMON /KONTRL/ DSCTOL, ERRORI, XLEFT, XRIGHT, BREAK, BOTH,
  1607.      * FACTOR, FMESGE, IBREAK, INTERP, LEFT, MAXAUX, MAXKNT, MAXPAR,
  1608.      * MAXSTK, NPAR, NSTACK, RIGHT, DISCRD, FATAL, FINISH
  1609.       COMMON /COMDIF/ DDTEMP, FINTRP, FLEFT, FRIGHT, XDD, XINTRP,
  1610.      * LEFTX, NINTRP, RIGHTX
  1611. C
  1612. C         SET SOME SHORT LOCAL VARIABLE NAMES
  1613.       XL = XLEFT(NSTACK)
  1614.       XR = XRIGHT(NSTACK)
  1615.       NL = LEFTX
  1616.       NL1 = NL + 1
  1617.       NR = RIGHTX
  1618.       NI = NINTRP
  1619.       NRL = NR + NL
  1620.       NRI = NR + NI
  1621.       NRI1 = NRI - 1
  1622.       NRLI = NRL + NI
  1623. C
  1624. C         STARTING REPRESENTATION IS (ASSUMING XL = 0 )
  1625. C
  1626. C    D(1) +D(2)X +D(3)X**2 + --- +D(NL)X**(NL-1)
  1627. C   +(X**NL)*( D(NL+1)(+D(NL+2)(X-XR)**2 + --- +D(NL+NR)*(X-XR)**(NR-1)
  1628. C        *((X-XR)**NR)*(D(NL+NR+1) + D(NL+NR+2)*(X-XINTRP(1))
  1629. C                      +D(NL+NR+3)*(X-XINTRP(1))(X-XINTRP(2)) + ---))
  1630. C
  1631. C         STRATEGY IS TO FIRST CONVERT THE PART FROM THE INTERP. PTS.
  1632. C         TO POLY IN (X-XR).  THIS POLY THEN HAS ORIGIN SHIFTED TO XL.
  1633. C
  1634. C     THE CONVERSION OF THE INTERP PART IS DONE EXPLICITLY FOR DEGREE
  1635. C     TWO OR LESS AND DONE BY SYNTHETIC DIVISION FOR HIGHER DEGREES
  1636. C
  1637. C   D1 + D2(X-X1) +D3(X**2-(X1+X2)X +X1*X2)
  1638. C
  1639. C     THE RESULTING COEFFICIENTS ARE PUT IN THE ARRAY POWERS
  1640. C
  1641. C             DEBUG  DEBUG  DEBUG  DEBUG
  1642.       IF (LEVEL.EQ.4) WRITE (6,99999) (D(K,1),K=1,NPAR)
  1643.       IF (NI.EQ.0) GO TO 100
  1644. C             BUILD UP THE POLYNOMIAL FOR THE INTERPOLATION POINTS
  1645. C
  1646. C         USE SPECIAL FORMULAS FOR NI LESS THAN 3
  1647.       IF (NI.EQ.1) GO TO 10
  1648.       IF (NI.EQ.2) GO TO 20
  1649.       GO TO 30
  1650.    10 POWERS(1) = D(NRL+1,1)
  1651.       GO TO 80
  1652.    20 POWERS(1) = D(NRL+1,1) + (XR-XINTRP(1))*D(NRL+2,1)
  1653.       POWERS(2) = D(NRL+2,1)
  1654.       GO TO 80
  1655. C
  1656. C         CONVERSION BY REPEATED SYNTHETIC DIVISION
  1657.    30 NI1 = NI - 1
  1658. C         INITIALIZE THE POWERS AND XTEMP ARRAYS
  1659.       DO 40 K=1,NI
  1660.         XTEMP(K) = XINTRP(K)
  1661.         NRLK = NRL + K
  1662.         POWERS(K) = D(NRLK,1)
  1663.    40 CONTINUE
  1664. C
  1665. C         DO THE REPEATED SYNTHETIC DIVISION TO REPLACE THE XTEMP
  1666. C         = XINTRP POINTS OF THE NEWTON EXPANSION BY THE XR POINTS
  1667.       DO 70 K=1,NI1
  1668. C              POWERS(NI) IS FIXED AND SET ABOVE
  1669.         DO 50 II=1,NI1
  1670.           I = NI - II
  1671.           POWERS(I) = POWERS(I) + (XR-XTEMP(I))*POWERS(I+1)
  1672.    50   CONTINUE
  1673. C              SHIFT THE NEWTON EXPANSION PTS. UP, PUT IN ONE MORE XR
  1674.         DO 60 II=1,NI1
  1675.           I = NI - II
  1676.           XTEMP(I+1) = XTEMP(I)
  1677.    60   CONTINUE
  1678.         XTEMP(1) = XR
  1679.    70 CONTINUE
  1680.       IF (LEVEL.EQ.4) WRITE (6,99998) (POWERS(K),K=1,NI)
  1681.    80 CONTINUE
  1682. C             SHIFT THE COEFFICIENTS TO THE TOP OF THE POWERS ARRAY
  1683.       DO 90 K=1,NI
  1684.         L = NI + 1 - K
  1685.         LTOP = L + NRL
  1686.         POWERS(LTOP) = POWERS(L)
  1687.    90 CONTINUE
  1688. C
  1689. C             HAVE THE INTERPOLATION PT. COEFS. IN THE ARRAY POWERS
  1690.   100 CONTINUE
  1691. C             PUT THE REMAINING DIVIDED DIFFS INTO THE POWERS ARRAY
  1692.       DO 110 J=1,NRL
  1693.         POWERS(J) = D(J,1)
  1694.   110 CONTINUE
  1695. C
  1696. C             DEBUG  DEBUG  DEBUG  DEBUG
  1697.       IF (LEVEL.EQ.4) WRITE (6,99997) (POWERS(K),K=NL1,NPAR)
  1698. C         TRANSFORM THE ORIGIN OF THE POLYNOMIAL FROM XR TO XL
  1699. C         WE USE REPEATED SYNTHETIC DIVISION
  1700.       IF (NRI.EQ.1) GO TO 140
  1701.       SHIFT = XR - XL
  1702.       KHI = NRI1
  1703. C             LOOP THROUGH THE COEFFICIENTS
  1704.       DO 130 J=2,NRI
  1705. C                   SYNTHETIC DIVISION LOOP
  1706.         DO 120 K=1,KHI
  1707.           KOEF = NRLI - K
  1708.           POWERS(KOEF) = POWERS(KOEF) - SHIFT*POWERS(KOEF+1)
  1709.   120   CONTINUE
  1710.         KHI = KHI - 1
  1711.   130 CONTINUE
  1712.   140 CONTINUE
  1713. C     THE COEFFICIENTS ARE NOW OF THE POWER FORM WITH ORIGIN XL
  1714. C             DEBUG  DEBUG  DEBUG  DEBUG
  1715.       IF (LEVEL.EQ.4) WRITE (6,99996) (POWERS(K),K=1,NPAR)
  1716.       RETURN
  1717. 99999 FORMAT (15X, 26HDEBUG PTRANS, ORIG INPUT =/(5X, 8E15.5))
  1718. 99998 FORMAT (10X, 17HINTERP PART COEFS/(5X, 8E15.5))
  1719. 99997 FORMAT (10X, 25HRIGHT + INTERP PART COEFS/(5X, 8E15.5))
  1720. 99996 FORMAT (10X, 18HFINAL COEFFICIENTS/(5X, 8E15.5))
  1721.       END
  1722.       SUBROUTINE PUT(INTERV, XKNOTS, COEFS, KDIMEN, NDIMEN)
  1723. C
  1724. C        ===============================================================
  1725. C **  THIS PROGRAM PUTS INTERVALS ON THE STACK OR DISCARDS THEM.
  1726. C     WHEN AN INTERVAL IS DISCARDED A NEW KNOT IS FOUND. THEN THIS
  1727. C     PROGRAM UPDATES THE ERROR ESTIMATE, THE XKNOT ARRAY, TRANSFORMS
  1728. C     THE POLYNOMIAL TO THE POWER FORM AND PUT THE COEFFICIENTS INTO
  1729. C     THE ARRAY COEFS.  IT ALSO CHECKS FOR PASSING BREAK POINTS
  1730. C
  1731.       REAL A, ACCUR, B, BLEFT, BRIGHT, CHARF, DDTEMP, DSCTOL, ERROR,
  1732.      * ERRORI, FACTOR, FINTRP, FLEFT, FRIGHT, NORM, XBREAK, XDD,
  1733.      * XINTRP, XLEFT, XRIGHT
  1734.       DIMENSION XBREAK(20), DBREAK(20), BLEFT(20), BRIGHT(20)
  1735.       DIMENSION XLEFT(50), XRIGHT(50), FACTOR(12), FMESGE(6)
  1736.       DIMENSION DDTEMP(12,12), FINTRP(10), FLEFT(6), FRIGHT(6),
  1737.      * XDD(12), XINTRP(10)
  1738.       INTEGER BOTH, BREAK, DBREAK, DEGREE, EDIST, FMESGE, RIGHT,
  1739.      * RIGHTX, SMOOTH
  1740.       LOGICAL DISCRD, FATAL, FINISH
  1741.       REAL XKNOTS(KDIMEN), COEFS(KDIMEN,NDIMEN)
  1742.       REAL BUFFER, DX, POWERS, P, RATIO
  1743.       DIMENSION POWERS(12)
  1744.       COMMON /INPUTZ/ A, B, ACCUR, NORM, CHARF, XBREAK, BLEFT, BRIGHT,
  1745.      * DBREAK, DEGREE, SMOOTH, LEVEL, EDIST, NBREAK, KNTDIM, NPARDM
  1746.       COMMON /RESULZ/ ERROR, KNOTS
  1747.       COMMON /KONTRL/ DSCTOL, ERRORI, XLEFT, XRIGHT, BREAK, BOTH,
  1748.      * FACTOR, FMESGE, IBREAK, INTERP, LEFT, MAXAUX, MAXKNT, MAXPAR,
  1749.      * MAXSTK, NPAR, NSTACK, RIGHT, DISCRD, FATAL, FINISH
  1750.       COMMON /COMDIF/ DDTEMP, FINTRP, FLEFT, FRIGHT, XDD, XINTRP,
  1751.      * LEFTX, NINTRP, RIGHTX
  1752. C
  1753.       DATA BUFFER /1.E-12/
  1754. C            CHECK FOR DISCARDING THE INTERVAL
  1755.       IF (DISCRD) GO TO 30
  1756. C         SUBDIVIDE INTERVAL AND PLACE ON STACK
  1757.       IF (NSTACK.LT.MAXSTK) GO TO 10
  1758. C             FATAL ERROR, EXCEEDED ACTIVE STACK SIZE
  1759.       FATAL = .TRUE.
  1760.       IF (LEVEL.GE.0) WRITE (6,99999) FMESGE, MAXSTK, XLEFT(NSTACK),
  1761.      * XRIGHT(NSTACK)
  1762.       DISCRD = .TRUE.
  1763.       GO TO 30
  1764. C
  1765.    10 DX = (XRIGHT(NSTACK)-XLEFT(NSTACK))*.5
  1766. C              CHECK FOR SMALL INTERVALS
  1767.       RATIO = DX/(ABS(A)+ABS(B))
  1768.       IF (RATIO.GT.BUFFER) GO TO 20
  1769.       DISCRD = .TRUE.
  1770.       FATAL = .TRUE.
  1771.       IF (LEVEL.GE.0) WRITE (6,99998) XLEFT(NSTACK), XRIGHT(NSTACK)
  1772.       GO TO 30
  1773.    20 CONTINUE
  1774.       NSTACK = NSTACK + 1
  1775.       XLEFT(NSTACK) = XLEFT(NSTACK-1)
  1776.       XLEFT(NSTACK-1) = XRIGHT(NSTACK-1) - DX
  1777.       XRIGHT(NSTACK) = XLEFT(NSTACK-1)
  1778. C             DEBUG  DEBUG  DEBUG  DEBUG
  1779.       IF (LEVEL.GE.3) WRITE (6,99997) NSTACK, XLEFT(NSTACK),
  1780.      * XRIGHT(NSTACK)
  1781.       RETURN
  1782. C
  1783. C            DISCARD INTERVAL, UPDATE GLOBAL ERROR, XKNOTS AND COEFS
  1784.    30 CONTINUE
  1785. C
  1786.       P = ABS(NORM)
  1787.       IF (NORM.EQ.3.) ERROR = AMAX1(ERROR,ERRORI)
  1788.       IF (NORM.NE.3.) ERROR = (ERROR**P+ERRORI)**(1./P)
  1789. C
  1790. C              CHECK FOR PASSING BREAK POINTS
  1791.       IF (BREAK.EQ.LEFT .OR. BREAK.EQ.BOTH) IBREAK = IBREAK + 1
  1792. C             DEBUG  DEBUG  DEBUG  DEBUG
  1793.       IF (LEVEL.GE.3) WRITE (6,99996) NSTACK
  1794. C
  1795. C             TRANSFORM REPRESENTATION OF POLYNOMIAL FROM DIVIDED
  1796. C             DIFFERENCES TO POWERS OF X WITH ORIGIN AT XKNOTS (KNOTS)
  1797.       CALL PTRANS(DDTEMP, POWERS)
  1798. C
  1799. C             PUT COEFS INTO THE MAIN ARRAY
  1800.       DO 40 K=1,NPAR
  1801.         COEFS(KNOTS,K) = POWERS(K)
  1802.    40 CONTINUE
  1803. C            PUT THE NEW KNOTS IN XKNOTS
  1804.       KNOTS = KNOTS + 1
  1805.       XKNOTS(KNOTS) = XRIGHT(NSTACK)
  1806.       NSTACK = NSTACK - 1
  1807.       RETURN
  1808. 99999 FORMAT (//2X, 6A4, 36HINTERVAL DIVIDED TOO MUCH, EXCEEDED ,
  1809.      * 6HLIMIT , I3, 22H ON INTERVAL STACK AT /20X, 2E18.8/2X, 6HINTERV,
  1810.      * 38HAL DISCARDED AND COMPUTATION CONTINUED)
  1811. 99998 FORMAT (25X, 23HGOT SHORT INTERVAL ****, 2E22.13, 12H **** DISCAR,
  1812.      * 4HD IT)
  1813. 99997 FORMAT (15X, 13HPUT INTERVAL , I3, 10H ON STACK , 2F12.8)
  1814. 99996 FORMAT (15X, 16HDISCARD INTERVAL, I5)
  1815.       END
  1816.       SUBROUTINE SETUP(XKNOTS, COEFS, KDIMEN, NDIMEN)
  1817. C
  1818. C        ===============================================================
  1819. C
  1820. C **  THIS PROGRAM CHECKS THE INPUT DATA AND INITIALIZES THE COMPUTATION
  1821. C
  1822.       REAL A, ACCUR, B, BLEFT, BRIGHT, CHARF, DDTEMP, DSCTOL, ERROR,
  1823.      * ERRORI, FACTOR, FINTRP, FLEFT, FRIGHT, NORM, XBREAK, XDD,
  1824.      * XINTRP, XLEFT, XRIGHT
  1825.       DIMENSION XBREAK(20), DBREAK(20), BLEFT(20), BRIGHT(20)
  1826.       DIMENSION XLEFT(50), XRIGHT(50), FACTOR(12), FMESGE(6)
  1827.       DIMENSION DDTEMP(12,12), FINTRP(10), FLEFT(6), FRIGHT(6),
  1828.      * XDD(12), XINTRP(10)
  1829.       INTEGER BOTH, BREAK, DBREAK, DEGREE, EDIST, FMESGE, RIGHT,
  1830.      * RIGHTX, SMOOTH
  1831.       LOGICAL DISCRD, FATAL, FINISH
  1832.       REAL XKNOTS(KDIMEN), COEFS(KDIMEN,NDIMEN)
  1833.       INTEGER HMSGE(6), NMSGE(4,4)
  1834.       COMMON /INPUTZ/ A, B, ACCUR, NORM, CHARF, XBREAK, BLEFT, BRIGHT,
  1835.      * DBREAK, DEGREE, SMOOTH, LEVEL, EDIST, NBREAK, KNTDIM, NPARDM
  1836.       COMMON /RESULZ/ ERROR, KNOTS
  1837.       COMMON /KONTRL/ DSCTOL, ERRORI, XLEFT, XRIGHT, BREAK, BOTH,
  1838.      * FACTOR, FMESGE, IBREAK, INTERP, LEFT, MAXAUX, MAXKNT, MAXPAR,
  1839.      * MAXSTK, NPAR, NSTACK, RIGHT, DISCRD, FATAL, FINISH
  1840.       COMMON /COMDIF/ DDTEMP, FINTRP, FLEFT, FRIGHT, XDD, XINTRP,
  1841.      * LEFTX, NINTRP, RIGHTX
  1842. C
  1843. C              IKNOT IS A COUNTER USED IN THE OUTPUT APPROX. PPOLY
  1844.       COMMON /SAVEIT/ IKNOT
  1845.       DATA HMSGE(1), HMSGE(2), HMSGE(3), HMSGE(4), HMSGE(5), HMSGE(6) /
  1846.      * 4H ***,4H* FA,4HTAL ,4HERRO,4HR **,4H**  /
  1847.       DATA NMSGE(1,1), NMSGE(1,2), NMSGE(1,3), NMSGE(1,4) /4HLEAS,
  1848.      * 4HT DE,4HVIAT,4HIONS/
  1849.       DATA NMSGE(2,1), NMSGE(2,2), NMSGE(2,3), NMSGE(2,4) /4HLEAS,
  1850.      * 4HT SQ,4HUARE,4HS   /
  1851.       DATA NMSGE(3,1), NMSGE(3,2), NMSGE(3,3), NMSGE(3,4) /4HMINI,
  1852.      * 4HMAX ,4HNORM,4H    /
  1853.       DATA NMSGE(4,1), NMSGE(4,2), NMSGE(4,3), NMSGE(4,4) /4HGENE,
  1854.      * 4HRAL ,4HL-P ,4HNORM/
  1855.       DATA KLEFT, KRIGHT, KBOTH /4HLEFT,4HRITE,4HBOTH/
  1856. C
  1857. C         PUT DATA STATEMENT ITEMS INTO COMMON VARIABLES
  1858. C         VARIOUS COMPILERS ALLOW MORE EFFICIENT WAYS
  1859. C
  1860.       LEFT = KLEFT
  1861.       RIGHT = KRIGHT
  1862.       BOTH = KBOTH
  1863.       DO 10 K=1,6
  1864.         FMESGE(K) = HMSGE(K)
  1865.    10 CONTINUE
  1866. C -------- SET CURRENT VALUES OF LIMITS ON DIMENSIONS ------------------
  1867.       KNTDIM = KDIMEN
  1868.       NPARDM = NDIMEN
  1869.       MAXKNT = KNTDIM
  1870.       MAXSTK = 50
  1871.       MAXPAR = MIN0(12,NPARDM)
  1872.       MAXAUX = 20
  1873. C -------- CHECK INPUT DATA --------------------------------------------
  1874.       FATAL = .FALSE.
  1875.       IF ((LEVEL+1)*LEVEL*(LEVEL-1)*(LEVEL-2).EQ.0) GO TO 20
  1876. C         FIXED ERROR ON PRINT CONTROL LEVEL
  1877.       WRITE (6,99999) LEVEL
  1878. C         FOR DEBUGGING DO NOT CHANGE OUTPUT LEVEL
  1879. C         FOR PRODUCTION VERSION SET LEVEL = 0
  1880.    20 IF (A.LT.B) GO TO 30
  1881. C         FATAL ERROR IN INTERVAL (A,B)
  1882.       FATAL = .TRUE.
  1883.       WRITE (6,99998) FMESGE, A, B
  1884.    30 IF (ACCUR.GT.0.0) GO TO 40
  1885. C         FATAL ERROR, NEGATIVE ACCURACY
  1886.       FATAL = .TRUE.
  1887.       WRITE (6,99997) FMESGE, ACCUR
  1888.    40 IF (DEGREE.LT.MAXPAR) GO TO 50
  1889. C         FATAL ERROR, DEGREE EXCEEDS MAXIMUM ALLOWED VALUE
  1890.       IF (LEVEL.GE.0) WRITE (6,99996) FMESGE, DEGREE
  1891.       FATAL = .TRUE.
  1892.    50 IF (2*SMOOTH.LT.DEGREE) GO TO 60
  1893. C         FATAL ERROR, SMOOTH AND DEGREE INCOMPATIBLE
  1894.       FATAL = .TRUE.
  1895.       WRITE (6,99995) FMESGE, SMOOTH, DEGREE
  1896.    60 IF (CHARF.GT.(B-A)/FLOAT(MAXKNT)) GO TO 70
  1897. C         FATAL ERROR, CHARF IS TOO SMALL
  1898.       FATAL = .TRUE.
  1899.       WRITE (6,99994) FMESGE, CHARF
  1900.    70 IF ((NORM-1.)*(NORM-2.)*(NORM-3.).EQ.0.0 .OR. NORM.LT.0.0) GO TO
  1901.      * 80
  1902. C         FATAL ERROR, UNDEFINED NORM REQUESTED
  1903.       IF (LEVEL.GE.0) WRITE (6,99993) FMESGE, NORM
  1904.       FATAL = .TRUE.
  1905.    80 IF (NBREAK.GE.0 .AND. NBREAK.LE.MAXAUX) GO TO 90
  1906. C         FATAL ERROR IN NUMBER OF BREAK POINTS
  1907.       FATAL = .TRUE.
  1908.       WRITE (6,99992) FMESGE, NBREAK
  1909. C
  1910.    90 IF (NBREAK.EQ.0) GO TO 150
  1911. C         CHECK THE BREAK POINT DATA, MONOTONICITY AND DEGREE
  1912.       J = 1
  1913.       IF (XBREAK(1).LT.A .OR. XBREAK(NBREAK).GT.B) GO TO 110
  1914.       IF (NBREAK.EQ.1) GO TO 120
  1915.       DO 100 J=2,NBREAK
  1916.         IF (XBREAK(J-1).GE.XBREAK(J)) GO TO 110
  1917.   100 CONTINUE
  1918.       GO TO 120
  1919.   110 FATAL = .TRUE.
  1920. C         BREAK POINTS ARE NOT MONOTONIC
  1921.       WRITE (6,99991) FMESGE, XBREAK(J)
  1922.   120 LIMSM = (DEGREE-1)/2
  1923.       DO 130 J=1,NBREAK
  1924.         IF (DBREAK(J).LT.0 .OR. DBREAK(J).GT.LIMSM) GO TO 140
  1925.   130 CONTINUE
  1926.       GO TO 150
  1927. C         BAD VALUE IN DERIVATIVE BREAKS
  1928.   140 FATAL = .TRUE.
  1929.       WRITE (6,99990) FMESGE, DBREAK(J)
  1930. C
  1931.   150 IF (EDIST*(EDIST-1)*(EDIST-2).EQ.0) GO TO 160
  1932. C         FATAL ERROR, ILLEGAL ERROR DISTRIBUTION REQUESTED
  1933.       IF (LEVEL.GE.0) WRITE (6,99989) FMESGE, EDIST
  1934.       FATAL = .TRUE.
  1935.   160 CONTINUE
  1936. C
  1937. C -------- INITIALIZATION OF VARIABLES ---------------------------------
  1938. C
  1939. C             ACTIVE INTERVAL STACK
  1940.       NSTACK = 1
  1941.       XLEFT(1) = A
  1942.       XRIGHT(1) = B
  1943. C             TERMINATION AND ERROR VALUES
  1944.       FINISH = .FALSE.
  1945.       ERROR = 0.
  1946.       DSCTOL = ACCUR**ABS(NORM)
  1947.       IF (EDIST.EQ.0) DSCTOL = DSCTOL/(B-A)
  1948.       IF (NORM.EQ.3.) DSCTOL = ACCUR
  1949. C             MISCELLANEOUS VARIABLES AND POINTERS
  1950.       IBREAK = 1
  1951.       KNOTS = 1
  1952.       IKNOT = 1
  1953.       INTERP = DEGREE + 2 - 2*SMOOTH
  1954.       XKNOTS(1) = A
  1955.       NPAR = DEGREE + 1
  1956. C
  1957. C             COMPUTE ARRAY OF NPAR FACTORIALS, FACTOR(K)= K-1 FACTORIAL
  1958.       FACTOR(1) = 1.
  1959.       FACTOR(2) = 1.
  1960.       DO 170 K=3,NPAR
  1961.         FACTOR(K) = FLOAT(K-1)*FACTOR(K-1)
  1962.   170 CONTINUE
  1963. C
  1964. C -------- PRINT OUT OF PROBLEM TO BE SOLVED ---------------------------
  1965.       IF (LEVEL.LE.0) GO TO 180
  1966.       NMES = NORM
  1967.       IF (NORM.LT.0.0) NMES = 4
  1968.       WRITE (6,99988) A, B, DEGREE, SMOOTH, ACCUR, (NMSGE(NMES,J),J=1,4)
  1969.       WRITE (6,99987) CHARF, NORM
  1970.       IF (EDIST.EQ.0) WRITE (6,99986)
  1971.       IF (EDIST.EQ.1) WRITE (6,99985)
  1972.       IF (EDIST.EQ.2) WRITE (6,99984)
  1973.       IF (EDIST.EQ.-1) WRITE (6,99983)
  1974.       IF (NBREAK.EQ.0) GO TO 180
  1975.       WRITE (6,99982) NBREAK
  1976.       WRITE (6,99981) (XBREAK(J),DBREAK(J),BLEFT(J),BRIGHT(J),J=1,
  1977.      * NBREAK)
  1978.   180 CONTINUE
  1979.       RETURN
  1980. 99999 FORMAT (//2X, 29HILLEGAL OUTPUT LEVEL CONTROL , I2, 9H SET TO 0)
  1981. 99998 FORMAT (//2X, 6A4, 20H INCORRECT INTERVAL , 2E15.5)
  1982. 99997 FORMAT (//2X, 6A4, 28H ILLEGAL ACCURACY REQUESTED , E15.5)
  1983. 99996 FORMAT (//2X, 6A4, 6HDEGREE, I3, 29H EXCEEDS MAXIMUM ALLOWED VALU,
  1984.      * 1HE)
  1985. 99995 FORMAT (//2X, 6A4, 35H INCOMPATIBLE DEGREE AND SMOOTHNESS, 2I5)
  1986. 99994 FORMAT (//2X, 6A4, 40H CHARACTERISTIC OSCILLATION LENGTH OF F ,
  1987.      * E15.5, 14H IS TOO SMALL )
  1988. 99993 FORMAT (//2X, 6A4, 20H ERROR MEASURE NORM , E15.5, 11H NOT ALLOWE,
  1989.      * 1HD)
  1990. 99992 FORMAT (//2X, 6A4, 11H IN NUMBER , I4, 17H OF BREAK POINTS )
  1991. 99991 FORMAT (//2X, 6A4, 30H BREAK POINTS NOT IN ORDER AT , E15.5)
  1992. 99990 FORMAT (//2X, 6A4, 20H ILLEGAL DERIVATIVE , I3, 10H AT BREAK )
  1993. 99989 FORMAT (//2X, 6A4, 32HILLEGAL ERROR DISTRIBUTION TYPE , I3)
  1994. 99988 FORMAT (//2X, 45H++++++ PIECEWISE POLYNOMIAL APPROXIMATION ON ,
  1995.      * 9HINTERVAL , 2E15.4, 11H OF DEGREE , I2, 6H WITH , I2, 7H CONTIN,
  1996.      * 17HUOUS DERIVATIVES /10X, 23H ACCURACY REQUESTED IS , E15.4,
  1997.      * 13H MEASURED BY , 4A4)
  1998. 99987 FORMAT (10X, 43HOTHER INPUT/DEFAULT VARIABLES ARE FOSCIL = ,
  1999.      * E15.4, 10H, EMEAS = , E15.4)
  2000. 99986 FORMAT (15X, 9(2H--), 33H PROPORTIONAL ERROR DISTRIBUTION )
  2001. 99985 FORMAT (15X, 9(2H--), 36HAPPROXIMATE FIXED ERROR DISTRIBUTION)
  2002. 99984 FORMAT (15X, 9(2H--), 25H FIXED ERROR DISTRIBUTION)
  2003. 99983 FORMAT (15X, 9(2H--), 28HMODIFIED PROPORTIONAL ERROR , 8HDISTRIBU,
  2004.      * 4HTION)
  2005. 99982 FORMAT (I12, 36H BREAK POINTS SPECIFIED WITH DATA = , 9H LOCATION,
  2006.      * 32H, DERIVATIVE, LEFT+RIGHT VALUES )
  2007. 99981 FORMAT (5X, 2(E20.5, I4, E16.5, E15.5))
  2008.       END
  2009.       SUBROUTINE SUMARY(XKNOTS, COEFS, KDIMEN, NDIMEN)
  2010. C
  2011. C        ===============================================================
  2012. C
  2013. C **  THIS PROGRAM PRINTS OUT A SUMMARY OF RESULTS OF ADAPT
  2014. C
  2015.       REAL A, ACCUR, B, BLEFT, BRIGHT, CHARF, DDTEMP, DSCTOL, ERROR,
  2016.      * ERRORI, FACTOR, FINTRP, FLEFT, FRIGHT, NORM, XBREAK, XDD,
  2017.      * XINTRP, XLEFT, XRIGHT
  2018.       DIMENSION XBREAK(20), DBREAK(20), BLEFT(20), BRIGHT(20)
  2019.       DIMENSION XLEFT(50), XRIGHT(50), FACTOR(12), FMESGE(6)
  2020.       DIMENSION DDTEMP(12,12), FINTRP(10), FLEFT(6), FRIGHT(6),
  2021.      * XDD(12), XINTRP(10)
  2022.       INTEGER BOTH, BREAK, DBREAK, DEGREE, EDIST, FMESGE, RIGHT,
  2023.      * RIGHTX, SMOOTH
  2024.       LOGICAL DISCRD, FATAL, FINISH
  2025.       REAL XKNOTS(KDIMEN), COEFS(KDIMEN,NDIMEN)
  2026.       COMMON /INPUTZ/ A, B, ACCUR, NORM, CHARF, XBREAK, BLEFT, BRIGHT,
  2027.      * DBREAK, DEGREE, SMOOTH, LEVEL, EDIST, NBREAK, KNTDIM, NPARDM
  2028.       COMMON /RESULZ/ ERROR, KNOTS
  2029.       COMMON /KONTRL/ DSCTOL, ERRORI, XLEFT, XRIGHT, BREAK, BOTH,
  2030.      * FACTOR, FMESGE, IBREAK, INTERP, LEFT, MAXAUX, MAXKNT, MAXPAR,
  2031.      * MAXSTK, NPAR, NSTACK, RIGHT, DISCRD, FATAL, FINISH
  2032.       COMMON /COMDIF/ DDTEMP, FINTRP, FLEFT, FRIGHT, XDD, XINTRP,
  2033.      * LEFTX, NINTRP, RIGHTX
  2034. C
  2035.       IF (LEVEL.EQ.-1) RETURN
  2036. C                   LEVEL 0 OUTPUT
  2037.       WRITE (6,99999) DEGREE, SMOOTH, KNOTS, ERROR
  2038.       IF (LEVEL.EQ.0) RETURN
  2039. C
  2040. C                   LEVEL 1 OUTPUT - KNOTS AND COEFFICIENTS
  2041.       KNOTS1 = KNOTS - 1
  2042.       IF (KNOTS.GT.15) GO TO 20
  2043. C
  2044. C               FORMAT FOR SMALL NO. OF KNOTS
  2045.       WRITE (6,99998)
  2046.       DO 10 K=1,KNOTS1
  2047.         WRITE (6,99997) K, XKNOTS(K), COEFS(K,1)
  2048.         WRITE (6,99996) (COEFS(K,J),J=2,NPAR)
  2049.    10 CONTINUE
  2050.       GO TO 40
  2051. C
  2052. C               FORMAT FOR LOTS OF KNOTS
  2053.    20 CONTINUE
  2054.       WRITE (6,99995)
  2055.       DO 30 K=1,KNOTS1
  2056.         WRITE (6,99994) K, XKNOTS(K), (COEFS(K,J),J=1,NPAR)
  2057.    30 CONTINUE
  2058.    40 RETURN
  2059. 99999 FORMAT (///48H --- ADAPTIVE PIECEWISE POLYNOMIAL APPROXIMATION,
  2060.      * 11H OF DEGREE , I2, 6H WITH , I2, 12H CONTINUOUS , 10HDERIVATIVE,
  2061.      * 8HS NEEDED, I4, 17H KNOTS FOR ERROR , E10.4)
  2062. 99998 FORMAT (8X, 13HKNOT LOCATION, 13X, 21HX-POWER COEFFICIENTS ,
  2063.      * 27HRELATIVE TO KNOT LOCATIONS )
  2064. 99997 FORMAT (I10, 2E20.12)
  2065. 99996 FORMAT (30X, E20.12)
  2066. 99995 FORMAT (3X, 39HK   K-TH INTERIOR KNOT     POWERS OF X , 7HRELATIV,
  2067.      * 20HE TO KNOT LOCATIONS )
  2068. 99994 FORMAT (I5, E25.12, 4E22.12/(E30.12, 4E22.12))
  2069.       END
  2070.       SUBROUTINE TAKE(INTERV)
  2071. C
  2072. C        ===============================================================
  2073. C
  2074. C **  THIS PROGRAM TAKES AN ACTIVE INTERVAL OFF THE TOP OF THE STACK
  2075. C     IT ALSO DOES MOST OF THE WORK OF LOCATING AND HANDLING
  2076. C     BREAK POINTS
  2077. C
  2078.       REAL A, ACCUR, B, BLEFT, BRIGHT, CHARF, DDTEMP, DSCTOL, ERROR,
  2079.      * ERRORI, FACTOR, FINTRP, FLEFT, FRIGHT, NORM, XBREAK, XDD,
  2080.      * XINTRP, XLEFT, XRIGHT
  2081.       DIMENSION XBREAK(20), DBREAK(20), BLEFT(20), BRIGHT(20)
  2082.       DIMENSION XLEFT(50), XRIGHT(50), FACTOR(12), FMESGE(6)
  2083.       DIMENSION DDTEMP(12,12), FINTRP(10), FLEFT(6), FRIGHT(6),
  2084.      * XDD(12), XINTRP(10)
  2085.       INTEGER BOTH, BREAK, DBREAK, DEGREE, EDIST, FMESGE, RIGHT,
  2086.      * RIGHTX, SMOOTH
  2087.       LOGICAL DISCRD, FATAL, FINISH
  2088.       REAL DX, RATIO, BUFFER
  2089.       COMMON /INPUTZ/ A, B, ACCUR, NORM, CHARF, XBREAK, BLEFT, BRIGHT,
  2090.      * DBREAK, DEGREE, SMOOTH, LEVEL, EDIST, NBREAK, KNTDIM, NPARDM
  2091.       COMMON /RESULZ/ ERROR, KNOTS
  2092.       COMMON /KONTRL/ DSCTOL, ERRORI, XLEFT, XRIGHT, BREAK, BOTH,
  2093.      * FACTOR, FMESGE, IBREAK, INTERP, LEFT, MAXAUX, MAXKNT, MAXPAR,
  2094.      * MAXSTK, NPAR, NSTACK, RIGHT, DISCRD, FATAL, FINISH
  2095.       COMMON /COMDIF/ DDTEMP, FINTRP, FLEFT, FRIGHT, XDD, XINTRP,
  2096.      * LEFTX, NINTRP, RIGHTX
  2097. C
  2098.       DATA BUFFER /1.E-12/
  2099. C
  2100. C         CHECK FOR BREAK POINT
  2101.       BREAK = 0
  2102.       IF (NBREAK.EQ.0 .OR. IBREAK.GT.NBREAK) GO TO 20
  2103.       IF (XBREAK(IBREAK).GT.XRIGHT(NSTACK)) GO TO 20
  2104. C
  2105. C               SET CONTROL VARIABLE BREAK, CHECK FOR LOCATION
  2106.       IF (XBREAK(IBREAK).GT.XLEFT(NSTACK)) GO TO 10
  2107.       BREAK = LEFT
  2108.       IF (IBREAK.EQ.NBREAK) GO TO 20
  2109. C             CHECK FOR SECOND BREAK POINT IN THIS INTERVAL
  2110.       IF (XBREAK(IBREAK+1).GE.XRIGHT(NSTACK)) GO TO 20
  2111. C                 NEXT BREAK IS INSIDE INTERVAL, SPLIT TOP INTERVAL
  2112.       BREAK = BOTH
  2113. C                  CHECK EXCEEDING STACK LIMIT. IF SO, STOP
  2114.       IF (NSTACK.EQ.MAXSTK) GO TO 30
  2115. C                     DONT SPLIT VERY SMALL INTERVALS, STOP INSTEAD
  2116.       DX = XBREAK(IBREAK+1) - XLEFT(NSTACK)
  2117.       RATIO = DX/(ABS(A)+ABS(B))
  2118.       IF (RATIO.LE.BUFFER) GO TO 40
  2119.       NSTACK = NSTACK + 1
  2120.       XLEFT(NSTACK) = XLEFT(NSTACK-1)
  2121.       XRIGHT(NSTACK) = XBREAK(IBREAK+1)
  2122.       XLEFT(NSTACK-1) = XRIGHT(NSTACK)
  2123.       IF (LEVEL.GE.2) WRITE (6,99999) NSTACK, XLEFT(NSTACK),
  2124.      * XRIGHT(NSTACK)
  2125.       GO TO 20
  2126. C
  2127.    10 BREAK = RIGHT
  2128. C                 CHECK TO SEE IF BREAK IS ALREADY AT RIGHT END POINT
  2129.       IF (XBREAK(IBREAK).GE.XRIGHT(NSTACK)) GO TO 20
  2130. C                 THE BREAK IS INSIDE INTERVAL, SPLIT TOP INTERVAL
  2131. C                  CHECK EXCEEDING STACK LIMIT. IF SO, STOP
  2132.       IF (NSTACK.EQ.MAXSTK) GO TO 30
  2133. C                     DONT SPLIT VERY SMALL INTERVALS, STOP INSTEAD
  2134.       DX = XBREAK(IBREAK) - XLEFT(NSTACK)
  2135.       RATIO = DX/(ABS(A)+ABS(B))
  2136.       IF (RATIO.LE.BUFFER) GO TO 40
  2137.       NSTACK = NSTACK + 1
  2138.       XLEFT(NSTACK) = XLEFT(NSTACK-1)
  2139.       XRIGHT(NSTACK) = XBREAK(IBREAK)
  2140.       XLEFT(NSTACK-1) = XRIGHT(NSTACK)
  2141.       IF (LEVEL.GE.2) WRITE (6,99998) NSTACK, XLEFT(NSTACK),
  2142.      * XRIGHT(NSTACK)
  2143.    20 CONTINUE
  2144. C             DEBUG  DEBUG  DEBUG
  2145.       IF (LEVEL.GE.3) WRITE (6,99997) XLEFT(NSTACK), XRIGHT(NSTACK),
  2146.      * BREAK
  2147.       RETURN
  2148. C
  2149. C        UNABLE TO PROCEED BECAUSE THE STACK LIMIT IS REACHED WITH
  2150. C        MULTIPLE BREAKPOINTS INSIDE THE SMALLEST ONE.
  2151.    30 FATAL = .TRUE.
  2152.       FINISH = .TRUE.
  2153.       IF (LEVEL.GE.0) WRITE (6,99996) FMESGE, MAXSTK, XLEFT(NSTACK),
  2154.      * XRIGHT(NSTACK), XBREAK(IBREAK)
  2155.       RETURN
  2156. C
  2157. C         SPLITTING INTERVALS TO ACCOMODATE BREAK POINTS HAS LED TO
  2158. C         AN EXCESSIVELY SMALL INTERVAL
  2159.    40 FATAL = .TRUE.
  2160.       FINISH = .TRUE.
  2161.       IF (LEVEL.GE.0) WRITE (6,99995) FMESGE, XLEFT(NSTACK),
  2162.      * XBREAK(IBREAK), XRIGHT(NSTACK)
  2163.       RETURN
  2164. 99999 FORMAT (10X, 28H++++ SPLIT TOP INTERVAL, GET, I3, 2E18.8,
  2165.      * 17H FOR BREAK = BOTH)
  2166. 99998 FORMAT (10X, 28H++++ SPLIT TOP INTERVAL, GET, I3, 2E18.8,
  2167.      * 19H FOR BREAK = RIGHT )
  2168. 99997 FORMAT (15X, 14HTAKE INTERVAL , 2F12.8, 11H    BREAK =, A4)
  2169. 99996 FORMAT (//2X, 6A4, 42H INTERVAL DIVIDED TOO MUCH, EXCEEDED LIMIT,
  2170.      * I3, 21H ON INTERVAL STACK AT/20X, 2E18.8/2X, 12HBREAK POINT ,
  2171.      * E18.8, 48H PRESENT REQUIRES SPLITTING INTERVAL, STOP ADAPT)
  2172. 99995 FORMAT (//2X, 6A4, 43HBREAK POINT ADJUSTMENT HAS LED TO A FATALLY,
  2173.      * 1H , 40HSMALL INTERVAL, THE POINTS INVOLVED ARE /20X, 3E18.8)
  2174.       END
  2175.       SUBROUTINE TERMIN(TEST, AND, PRINT, XKNOTS, KDIMEN)
  2176. C
  2177. C        ===============================================================
  2178. C
  2179. C **  THIS PROGRAM TESTS FOR TERMINATION AND GIVES INTERMEDIATE OUTPUT
  2180. C
  2181.       REAL A, ACCUR, B, BLEFT, BRIGHT, CHARF, DDTEMP, DSCTOL, ERROR,
  2182.      * ERRORI, FACTOR, FINTRP, FLEFT, FRIGHT, NORM, XBREAK, XDD,
  2183.      * XINTRP, XLEFT, XRIGHT
  2184.       DIMENSION XBREAK(20), DBREAK(20), BLEFT(20), BRIGHT(20)
  2185.       DIMENSION XLEFT(50), XRIGHT(50), FACTOR(12), FMESGE(6)
  2186.       DIMENSION DDTEMP(12,12), FINTRP(10), FLEFT(6), FRIGHT(6),
  2187.      * XDD(12), XINTRP(10)
  2188.       INTEGER BOTH, BREAK, DBREAK, DEGREE, EDIST, FMESGE, RIGHT,
  2189.      * RIGHTX, SMOOTH
  2190.       LOGICAL DISCRD, FATAL, FINISH
  2191.       REAL XKNOTS(KDIMEN)
  2192.       COMMON /INPUTZ/ A, B, ACCUR, NORM, CHARF, XBREAK, BLEFT, BRIGHT,
  2193.      * DBREAK, DEGREE, SMOOTH, LEVEL, EDIST, NBREAK, KNTDIM, NPARDM
  2194.       COMMON /RESULZ/ ERROR, KNOTS
  2195.       COMMON /KONTRL/ DSCTOL, ERRORI, XLEFT, XRIGHT, BREAK, BOTH,
  2196.      * FACTOR, FMESGE, IBREAK, INTERP, LEFT, MAXAUX, MAXKNT, MAXPAR,
  2197.      * MAXSTK, NPAR, NSTACK, RIGHT, DISCRD, FATAL, FINISH
  2198.       COMMON /COMDIF/ DDTEMP, FINTRP, FLEFT, FRIGHT, XDD, XINTRP,
  2199.      * LEFTX, NINTRP, RIGHTX
  2200. C
  2201. C         INTERMEDIATE OUTPUT - ACCORDING TO LEVELS
  2202. C         NO INTERMEDIATE OUTPUT FOR LEVELS -1,0,1
  2203.       IF (LEVEL.LE.1) GO TO 40
  2204. C
  2205. C             LEVEL 2 OUTPUT - ONLY FOR DISCARDED INTERVAL
  2206.       IF (.NOT.DISCRD) GO TO 10
  2207.       WRITE (6,99999) KNOTS, XKNOTS(KNOTS), ERRORI, ERROR
  2208.       IF (BREAK.EQ.LEFT) WRITE (6,99998)
  2209.       IF (BREAK.EQ.RIGHT) WRITE (6,99997)
  2210.       IF (BREAK.EQ.BOTH) WRITE (6,99996)
  2211. C
  2212.       GO TO 20
  2213. C             DEBUG OUTPUT  -  LEVEL 3
  2214.    10 IF (LEVEL.EQ.2) GO TO 40
  2215. C
  2216. C                  INTERVAL SUMMARY
  2217.       WRITE (6,99995) NSTACK, XLEFT(NSTACK), XRIGHT(NSTACK), ERRORI
  2218.       IF (BREAK.NE.0) WRITE (6,99994) IBREAK, XBREAK(IBREAK)
  2219. C
  2220. C             DEBUG POLYNOMIAL OPERATIONS  -  LEVEL 4
  2221.    20 CONTINUE
  2222.       IF (LEVEL.LE.3) GO TO 40
  2223.       WRITE (6,99993) LEFTX, RIGHTX, NINTRP
  2224.       DO 30 K=1,NPAR
  2225.         KNPAR = NPAR + 1 - K
  2226.         WRITE (6,99992) K, (DDTEMP(I,K),I=1,KNPAR)
  2227.    30 CONTINUE
  2228.    40 CONTINUE
  2229. C             TEST FOR NORMAL TERMINATION
  2230.       IF (NSTACK.EQ.0) GO TO 50
  2231. C
  2232. C             TEST FOR ABNORMAL TERMINATION
  2233.       IF (KNOTS.LT.MAXKNT) RETURN
  2234. C
  2235. C            HAVE EXCEEDED LIMIT ON KNOTS
  2236.       WRITE (6,99991) FMESGE, MAXKNT, XKNOTS(KNOTS), ERROR
  2237.       FATAL = .TRUE.
  2238. C
  2239. C         TERMINATE COMPUTATION
  2240.    50 FINISH = .TRUE.
  2241.       RETURN
  2242. 99999 FORMAT (2X, 9H**** KNOT, I4, 4H AT , E16.5, 17H, WITH LOCAL AND ,
  2243.      * 16HGLOBAL ERRORS = , 2E15.4)
  2244. 99998 FORMAT (20X, 32HBREAK POINT ON LEFT OF INTERVAL )
  2245. 99997 FORMAT (20X, 33HBREAK POINT ON RIGHT OF INTERVAL )
  2246. 99996 FORMAT (20X, 37HBREAK POINT AT BOTH ENDS OF INTERVAL )
  2247. 99995 FORMAT (10X, 19H++++ STACK ELEMENT , I3, 3H = , E20.5, E15.5,
  2248.      * 19H, WITH LOCAL ERROR , E15.4, 17H PLACED ON STACK )
  2249. 99994 FORMAT (15X, 14HBBBREAK POINT , I3, 3H = , E15.4, 9H INVOLVED)
  2250. 99993 FORMAT (5X, 44H---- DIVIDED DIFFERENCE INFO ---- NL,NR,NI =, 3I3)
  2251. 99992 FORMAT (5X, 12HDIV DIFF ROW, I3, 9F12.5/93X, 3F12.5)
  2252. 99991 FORMAT (//2X, 6A4, 14HEXCEEDED LIMIT, I3, 14H ON NUMBER OF ,
  2253.      * 9HKNOTS AT , E18.8, 13H WITH ERROR =, E14.4)
  2254.       END
  2255.  
  2256.